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ABSTRACT 


In  this  network  synthesis  problem  a  matrix  giving  flow 
requirements  between  each  pair  of  points  is  specified,  and 
the  cost  of  flow  in  each  arc  is  a  concave  function  of  the 
amount  of  flow.  A  flow  pattern  which  fulfills  the  require¬ 
ments  at  minimum  cost  is  sought.  This  problem  is  formulated 
as  a  concave  programming  problem  with  linear  constraints. 

All  the  practical  difficulties  of  formulation  and  theoretical 
difficulties  of  identifying  the  globally  minimal  solution 
while  avoiding  locally  minimal  solutions  are  discussed. 

Three  procedures  are  developed.  The  first  procedure  is 
inefficient.  The  second  procedure  is  quite  efficient  but 
produces  locally  minimal  solutions.  Therefore  different 
heuristics  and  sequential  sampling  plans  are  suggested  to 
obtain  the  globally  minimal  solution  (or  estimate  it)  from  a 
number  of  locally  minimal  solutions.  The  third  procedure  is 
a  global  search  procedure  of  branch  and  bound  type  where  each 
subproblem  is  easy  to  solve.  This  procedure  has  limitations 
because  computer  core  size  requirement  grows  very  rapidly 
with  the  size  of  the  network.  Different  ways  of 
circumventing  this  limitation  are  discussed.  Finally,  a 
post  optimization  procedure  that  approximates  the  solution 
of  nonconcave  problems  with  cost  specified  by  step  functions 
defined  only  at  integral  points  is  discussed.  In  conclusion, 
a  more  general  problem  and  future  research  directions  are 
mentioned. 


CHAPTER  1 


1.1  Introduction 

This  thesis  studies  the  techniques  of  planning  a  communication 
network.  It  is  a  synthesis  problem;  the  required  number  of  channels 
connecting  each  source-sink  pair  is  specified  and  the  layout  of  the 
channels  of  communication  must  be  determined.  The  problem  is  sometimes 
known  as  the  "TELPAK"  problem  which  refer;;  to  packing  more  channels 
together  to  take  advantage  of  economies  of  scale. 

This  synthesis  problem  arises  most  often  in  two  situations.  The  first 
is  in  planning  an  intercity  network  such  that  the  cost  incurred  is 
minimum  given  an  estimate  of  the  traffic  between  cities.  The  traffic 
volume  data  are  used  in  C.C.S.^  calculations  to  determine  the  number  of 
channels  required  between  cities  for  a  certain  prespecified  probability  of 
blocking  and  certain  assumptions  about  the  alternate  routing  procedure. 

The  second  is  in  determining  the  most  economical  way  of  leasing,  from  a 
telephone  company,  trunk  groups  between  the  offices  of  a  corporation  given 
an  estimate  of  the  traffic  volume  between  the  offices  and  taking  advantage 
of  high  volume  discounts. 

In  both  the  above  cases,  the  planning  problem  would  be  fairly  simple, 
as  we  explain  later  in  this  chapter,  if  the  cost  per  channel  were  constant 
or  increasing  with  the  number  of  channels.  But  in  this  case,  we  face  a 
situation  where  improved  technology  can  be  used  when  the  number  of  channels 
increases.  Cost  per  channel  decreases  with  increasing  capacity.  This 
trend  is  evident  in  the  past  developments  (open  wire,  N2  carrier,  microwave 
radre,  analog  transmission  over  coaxial  cable)  and  is  expected  to  continue 
in  the  foreseeable  future  (waveguides). 

^C.C.S.  means  hundred  call  second  which  is  a  unit  used  in  measuring  traffic 
in  a  telephone  network. 
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When  the  cost  per  channel  is  a  decreasing  function,  we  have  a 
"concave  type"  (explained  in  the  note  on  Page  8)  function  as  the  total 
cost  objective  function  which  has  ;o  be  minimized.  A  small  example  of 
such  a  problem  should  be  useful.  Assume  that  a  single  telephone  line 
costs  one  dollar  per  mile,  but  if  we  build  one  hundred  lines  between  two 
points,  the  cost  is  75  dollars  per  mile.  If  we  build  two  hundred  lines 
between  two  cities,  it  costs,  say,  120  dollars  per  mile.  Because  the 
cost  is  concave  type,  it  is  cheaper  to  pack  the  lines  into  major  routes 
than  to  build  all  telephone  lines  directly.  To  illustrate  the  problem, 
suppose  that  we  have  six  locations  A  ,  B  ,  C  ,  D  ,  E  ,  and  F  as  shown 
in  Figure  2.  The  numbers  in  Figure  2  are  the  mileages  between  tue  points. 
Let  there  be  50  telephone  lines  required  between  points  A  and  C  ,  and 
also  between  points  B  and  D  .  No  lines  are  required  between  any  other 
pair  of  points.  The  total  cost  of  building  a  direct  line  between  A  and 
C  and  also  between  B  and  D  is 

2  x  (20  x  1  x  50)  =  2,000  dollars. 

If  we  do  not  build  direct  lines  between  B  and  D  but  build  the 
50  lines  by  way  of  A  and  C  ,  then  the  total  number  of  lines  to  be 
built  between  A  and  C  is  100,  and  the  cheaper  rate  is  available.  Thus, 
the  total  cost  is 

2  x  4  x  1  x  50  +  100  x  .75  x  20  =  1,900  dollars. 

If  we  do  not  build  direct  lines  between  either  A  and  C  or 
between  B  and  D  but  build  50  lines  between  A  and  E  ,  B  and  E  , 

C  and  F  ,  D  and  F  ,  and  also  100  lines  between  E  and  F  ,  then  the 


total  cost  is 
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4  x  2.8  x  1  x  50  +  100  x  .75  x  16  =  1,760  dollars. 

So,  the  natural  question  is:  what  is  the  cheapest  network  that 
satisfies  all  the.  requirements.  In  the  above  example,  the  third  type  is 
the  cheapest  so  far  obtained  and  to  get  that  we  had  to  enumerate  the 
cost  of  the  direct  routing  and  the  cost  of  using  other  possible  inter¬ 
mediate  nodes.  This  enumeration  becomes  difficult  as  the  number  of 
sources  and  sinks  increases  or  as  the  total  number  of  nodes  increases, 
and  becomes  almost  impossible  when  the  network  conf iguration  gets  more 
and  more  complex.  An  attempL  is  made  in  this  thesis  to  develop  a 
systematic  procedure  for  solving  moderate-sized  problems  of  this  type. 

1.2  Difficulties  Faced  in  Solving  the  P roblem 

The  difficulties  encountered  in  this  problem  can  be  roughly  divided 
into  two  categories.  The  first  type  of  difficulty  arises  in  building 
(from  a  real  world  communication  network  situation)  a  network  model  which 
can  be  analyzed.  This  will  be  illustrated  with  the  network  description 
that  follows  but  not  discussed  further  in  this  thesis.  The  second  type 
of  difficulty  arises  in  the  analysis  of  this  problem.  Because  of  a 
certain  mathematical  subtlety  involved  in  this  problem,  a  theoretically 
elegant  solution  procedure  appears  difficult  to  find,  if  one  exists  at 
all.  This  mathematical  difficulty  will  he  described  at  the  end  of  this 
section  to  justify  the  use  of  a  heuristic  solution  method  in  Chapter  3 
and  an  elaborate  enumerative  solution  method  in  Chapter  4. 

1  ■  ?.  ■  1  Network  Description  and_  Practical  Difficulties 

Telephone  calls  are  generated  by  the  individual  subscriber.  The 
subscriber's  demands  on  the  network  can  occur  whenever  he  choose;;,  can  he 
of  any  duration,  and  can  be  directed  to  any  other  subscriber.  Calls  may 
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be  rcuted  on  one  or  more  of  three  overlapping  networks.  The  network 
cc-nccting  each  subscriber  to  a  central  switching  of' ice  is  called  the 
local  network.  The  network  connecting  central  offices  within  an  area  of 
one  or  more  states  is  called  the  statewide  network.  The  network  which 
connects  the  main  cities  is  called  the  long-haul  network.  These  three 
classes  of  networks  cover  the  country.  The  synthesis  methods  developed 
in  this  thesis  are  applicable  mainly  to  the  long-haul  network  and  state¬ 
wide  networks.  In  the  case  of  statewide  networks,  the  main  stations  are 
the  nodal  points  of  the  system.  In  the  case  of  the  long-haul  network, 
cities  are  the  nodal  points  of  the  system. 

To  calculate  the  number  of  channels  required  between  nodal  points,  a 
projection  is  first  made  about  traffic  demand  between  two  points.  It  is 
very  difficult  to  project  these  traffic  demands  ahead  of  time.  Demand 
can  be  affected  by  a  change  in  rate  structure!.,  new  service  offerings,  or 
conditions  imposed  by  regulatory  agencies.  Moreover,  new  and  unexpected 
technologies  may  significantly  change  the  demand,  as,  for  example,  a 
communication  satellite  which  could  stimulate  long  distance  calling  and 
on-the-spot  T.V.  coverage,  or  an  improved  time  sharing  computer  system 
which  could  increase  data  transmission  through  communication  cables.  How¬ 
ever,  despite  all  these  uncertainties,  projections  are  made  for  traffic 
requirements  and  these  uncertainties  arc  taken  into  account.  Once  traffic 
projections  are  made,  the  C.C.S.  calculations  arc-  used  to  find  the  number 
of  channels  required  for  a  certain  probability  of  blocking.  Usually,  the 
blocking  probability  is  .01  and  standard  traffic  engineering  practices  are 
available  for  such  calculations. 

At  this  point,  two  complexities  whi<h  enter  into  trunk  requirement 
calculations  must  be  taken  into  account.  The  first  is  the  alternate 


routing  plan. 


The  alternate  route  is  the  pat  i  offered  to  a  call  if  the 


most  direct  path  is  blocked.  Alternate  routes  are  always  specified.  If 
an  alternate  route  from  point  A  to  point  B  goes  through  point  C  , 
then  this  fact  must  be  taken  into  account  in  determining  trunk  require¬ 
ments  between  points  A  and  C  ,  and  between  points  B  and  C  .  The 
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second  complexity  is  what  is  known  as  diversity  routing.  The  trunks 
between  two  points  are  split  into  two  or  more  physically  different  path? 
for  better  transmission  reliability.  That  is,  two  points  should  not  be 
completely  disconnected  if  a  part  of  the  network  is  destroyed  by  a  natural 
catastrophe  or  for  any  other  reason.  This  is  similar  to  alternate  routing 
requirements.  The  above  discussions  show  how  difficult  it  is  for  a  planner 
to  say  exactly  how  nany  channels  are  needed  between  each  pair  of  points. 
However,  for  our  synthesis  problem,  we  shall  assume  that  this  figure  has 
been  determined. 

Cost  Estimation 

There  is,  however,  another  set  of  difficulties,  even  given  the 
knowledge  of  these  channel  requirements.  Each  traffic  generating  and 
terminating  point,  and  all  traffic  handling  points,  are  defined  as  nodes, 
and  the  links  or  connecting  arcs  between  two  nodes  are  transmission 
facilities.  The  graph  is  undirected  because  telephone  conversations  go 
both  ways.  Several  facilities  may  comprise  a  link.  This  type  of  mapping 
by  means  of  nodes  and  arcs  is  not  always  simple — a  lot  of  judgemental 
factors  are  involved  in  deciding  these  nodes  and  arcs.  For  instance,  at 
certain  nodes  some  transmission  lines  pass  directly  through  and  some  go 
via  a  switching  mach*ne.  In  this  situation,  the  cost  is  different  for  the 
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See  "The  Design  of  Minimum  Cost  Survival  Networks,"  K.  Steiglitz, 

P.  Weiner,  D.  J.  Kleitman,  Technical  Memorandum  TM-105,  National  Resource 
Analysis  Center,  System  Evaluation  Division.  This  paper  gives  a 
heuristic  treatment. 


v  ,  c  to  the  extra  erst  of  switching.  At  other  nodes,  ail  calls 

ho  cw  tc  ''■•I.  7n  ?nv  case,  however,  the  switching  cost  can  be  taken 
•.  ■  _rr.0 Ky  a  technique  called  node  splitting  which  weans  that 

o-c  traffic  hand’ - ng  point  is  split  into  two  nodes  and  the  cost  of  the 
arc  between  .her-  is  'he  switching  cost.  So,  with  imagination  some  of  the 
d.‘  ff  icui  Lies  can  he  taken  care  of.  But  a  large  number  of  practical 
difficulties  cannot  be  tackled.  Experience  and  common  sense  must  be  used 
in  such  situations. 

Last,  but  not  the  least,  of  the  first  category  of  difficulties  is 
determining  the  cost  functions  which  relate  cost  to  the  number  of  channels. 
Cost  is  the  main  objective  in  our  design  considerations.  But  cost  to 
whom,  company  or  subscriber?  And  how  much  is  the  cost  incurred?  These 
are  natural  questions  which  arise  when  we  are  talking  about  cost.  In  the 
TELPAK  problem  where  interoffice  facilities  for  a  large  corporation  are 
considered,  the  cost  referred  to  is  cost  to  the  customer,  and  it  could  he 
different  from  the  actual  cost  incurred  by  the  telephone  company.  But  in 
other  cases  cost  refers  to  the  cost  incurred  by  the  telephone  company.  In 
both  situations,  an  element  of  uncertainty  is  involved  in  these  costs.  It 
is  affected  by  the  interest  rate,  available  funds  and  the  production  cost 
of  facilities  and  designs. 

The  possible  kinds  of  cost  functions  encountered  in  practice  in  this 

problem  are  illustrated  in  Figures  1.3,  1.4,  1.5,  1.6,  and  they  all  have 

one  thing  in  common — economy  of  scale.  All  the  cost  functions  are  what  we 
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will  call  "concave  type"  and  they  are  nondecreasing  functionsof  flow  value. 
_ 

That  is,  they  are  not  necessarily  concave  because  concave  f'inctions  are 
continuous  everywhere  except  at  the  end  points  and  the  following  functions 
may  not  have  that  continuity  property.  Also,  concave  functions  satisfy 
A  * r  (x  j  )  (1  -  A ) • f (x^)  '  f  (Ax^  +  (1  -  A ) x ^ )  ,  0  A  1  ,  which  only  the 

third  and  fourth  of  our  functions  satisfies  (assuming  A  be  such  that 
(Ax^  +  (l  -  A)x,,)  is  a  point  at  which  I  i;  defined). 


I 


A  j 


Number  of  Channels  — *- 

FIGURE  1.3 


The  step  function  with  decreasing  step  size — this  type  of  cost 
structure  may  be  encountered  because  of  the  uses  of  different  types  of 
facilities  for  different  ranges  of  the  number  of  channels.  A  more 
sophisticated  carrier  can  be  used  when  the  number  of  channels  is  large, 
resulting  in  a  lesser  cost  per  channel  and  hence  smaller  step  sizes.  The 
change  of  type  of  carrier  ran  only  occur  at:  certain  volumes  of  channel 
requirements  resulting  in  jumps  occurring  in  the  cost  curve. 


Number  of  Channels  — 

FI  Cl  IRK  1  .  A 

Semi-step  f  mic  t  i  on-  - 1  In  steps  are  connected  by  slant'd  line?;.  The 
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steps  of  decreasing  size  occur  due  to  the  same  reason  as  above.  Here  a 
variable  cost  is  incurred  in  transitional  capacity  levels  of  different 
facilities,  i.e.,  a  facility  can  be  overloaded  at  some  additional  cost  per 
channel. 


Number  of  Channels  — ► 


FIGURE  1.5 


Piecewise  linear  concave  function — here  each  part  is  linear  but  the 
slope  decreases  with  increasing  flow  value.  Here  a  variable  cost  is 
incurred  for  each  additional  channel  but  this  variable  cost  depends  on  the 
range  it  is  in. 


Number  of  Channels  ■  ► 


FIGURE  1.6 

The  final  type  has  decreasing  forward  di I  Terences ,  i.e.,  if  f  (n)  Is 
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the  total  cost  of  n  channels  then  f (n  +  3 )  -  f (n)  >  f (n  +  2)  -  f (n  +  1)  > 
...  >  0  .  Here  total  cost  increases  with  the  addition  of  each  channel, 
but  cost  per  channel  decreases.  This  is  strictly  a  concave  function  of 
discrete  type.  In  cases  where  only  approximate  cost  data  is  available  and 
some  guess  work  is  involved  in  getting  the  cost  values,  such  functions  arc- 
used. 

This  concludes  our  discussion  of  the  practical  difficulties  which 
arise  in  the  problem  formulation  stage.  V.'e  assume  that  these  difficulties 
can  be  surmounted  and  that  the  cost  function,  like  the  channel  require¬ 
ments,  is  available  to  the  network  synthesizer. 

Discrete  vs.  Continuous  Values 

The  cost  data  discussed  above  is  discrete,  i.e.,  values  art  only 
available  for  integral  points.  For  the  purpose  of  methods  used  in 
Chapters  3  and  4,  we  will  employ  a  continuous  concave  representation  of 
the  data,  obtained  by  using  a  continuous  curve  through  the  discrete- 
points.  However,  such  concave  curves  cannot  be  drawn  through  the  discrete 
points  for  Figures  J.3  and  1.4.  Hence,  some  approximation  is  necessary. 

Use  of  such  a  continuous  concave  version  is  justified  (except  Figures  1.3 
and  1.4)  because  the  answer  we  get  is  always  integral  (as  shown  in 
Appendix  C) ,  and  the  values  of  the  data  at  those  points  match  the 
original  data.  For  Figures  1.3  and  1.4,  some  post  opt  imization  procedure 
is  needed.  Fxcept  in  Chapter  T>,  we  will  always  be  com  clued  with 
continuous  variables  and  concave  functions.  Another  a  1  (  o !  u.it  i  vo  way  of 
joining  the  discrete  points  is  by  piecewise  linear  cuivcs.  The  ulgoi i t hm 
described  in  Chapter  3  is  applicable  to  .such  curves  with  only  inim-i 
modifications,  as  mentioned  in  that  chapter. 
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1.2.2  Theoretical  Difficulties 


We  are  now  in  a  position  to  discuss  the  second  type  of  difficulties, 
those  which  arise  in  the  analysis  of  the  problem.  Here  we  will  consider  a 


network  with  N  nodes  and  M  arcs  (M  can  be  different  from 


N(N  -  1) 


because  some  ~s  between  nodes  are  not  possible)  joining  the  nodes 

through  which  calls  are  channeled.  The  requirement  matrix  R  =  {r^}  is 

specified  where  r  =  number  of  channels  required  between  node  i  and 

i  .  Also  the  function  f  (y  )  is  given  for  all  M  arcs.  This  function 

gives  the  cost  of  y  channels  (amount  of  flow,  in  network  terminology) 

"  m 

in  arc  m  of  the  network.  With  all  this  information,  we  can  set  up  the 
above  problem  as  a  programming  problem. 

In  the  discussion  that  follows,  the  concept  of  extreme  points  of  a  set 
is  very  important. 


Definitions: 

A  convex  combination  of  two  points  and  in  a  vector  space  is 

a  point  X  given  by  X  =  (XX^  +  (1  -  X)X£)  where  X  is  a  scalar  and 
0  <  A  <  1  . 

An  extreme  point  X  of  a  set  is  a  point  in  the  set  such  that  there 
exist  no  two  distinct  points  X^  and  X2  in  the  set  and  1  >  X  >  0  such 
that  the  point  X  is  a  convex  combination  of  X^  and  X2  • 

Extreme  flow  is  similarly  defined  as  a  feasible  flow  which  cannot  be 
expressed  as  a  convex  combination  of  two  other  distinctly  different 
feasible  flows. 

In  the  arc-chain  formulation  that  follows,  two  types  of  vector  spaces, 
X  and  Y  ,  are  used.  The  X-space  is  defined  as  "chain  spa''’"  and  each  of 
its  element  |xij|  indicate  the  amount  of  flow  in  the  kth  chain  between 
source  and  sink  pair  i  j  .  The  Y-space  is  defined  as  "arc-space"  and 


J 


1  3 


each  of  its  elements  {y  }  indicates  the  amount  of  flow  in  arc  m  . 

m 

Y-space  is  a  linear  transform  of  X-space  by  means  of  arc-chain  incidence 
matrices  defined  below. 


Arc-Chain  Formulation 


=  Number  of  different  loopless  chains  between  node  i  and  j  . 
=  The  set  of  loopless  chains  connecting  node  i  and  j  ,  where 


V  s  1  2  P 

f  >  •  •  •  t  * 


X. .  =  A  vector  of  dimension  P..  where  element 
ij  ij 


R} 


represents 
k 


x^  amount  of  flow  between  i  and  j  through  chain  p  . 

A..  =  M  *  P..  arc-chain  incidence  matrix,  whose  (m,k)  element  is 
ij  ij 

k 

1  if  the  chain  p..  traverses  arc  m  ,  or  0  if  it  does  not. 

Y  =  y  A..X..  =  M  vector  {y  }  which  is  the  total  flow  ir  the 
.  li  xi  m 

arc  m  . 


So,  here  we  have  to  minimize  the  cost  Z  ,  i.e.. 


(1.1a) 


M 

Minimize  Z  =  )  f  (y  )  , 

S  m  m 
m=l 


subject  to  the  constraints 


ii  k 

(1.1b)  )  x..  >  r..  for  all  source  and  sink  pairs  i  j  , 

k=l  1J  ”  3J 


(J.Jt) 


Xj  j  >_  0  for  all  k  and  all  source  and  sink  pairs  i  j 


T!ie  linear  inequalities  (1.1b)  and  (1.1c)  define  a  polyhedron  In 
X-space.  Since  Y-space  is  a  linear  transform  of  X-spacc,  the  above 
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polyhedron  can  be  transformed  into  another  polyhedron  in  Y-space.  The 
objective  function  Z  is  the  positive  sum  of  concave  functions  of 
elements  of  the  vector  Y  .  So,  it  is  a  concave  function  in  Y- 

space. 


m 


Z(Y)  “  zi  l  A,  X  \  de{inJ- 
\i,j?i>j  lj  lj/ 

objective  function  in  X-space. 


tlon  Z(X)  ,  where  Z(X)  is  the 


Z(X)  is  a  concave  function  in  X-space  as  evident  from  the  following 
•  1  2 

inequality.  Consider  any  two  points  X  ,  X  ,  in  X-space  and  a  scalar 

12  1 
A  such  that  0  A  1  .  Y  and  Y  are  defined  as  the  images  of  X 

and  X2  in  Y-space,  i.e.  ,  Y1  =  £  A  X*  and  Y2  =  £  . 


ij  ij 


Z(AX  +  (1 


- x,x2)  ■  z(Uj L >s  auK +  °  -  x>xi3)) by 


above  definition 


because  A  is  a  scalar  :nd  A^'s  are  matrices. 


-  Z[A 

V 


l  AX*  +  (1  -  A) 


j*i>j 


ij  ij 


l  A  X2.\ 

i.jii>J  ij  1J/ 


=  Z(AY*  +  (1  -  A)Y2)  by  definition 

>_  AZ(Y*)  +  (1  -  A)Z(Y2)  because  Z(Y)  is  concave 
in  Y-space 


(i  Jiu  ^  +  a  '  X>i aU  V«)  by 


definition 


-  AZ(X*)  +  (1  -  A)*Z(X2)  by  definition 


The  following  lemma  is  useful  in  determining  extreme  points  of  the 
polyhedral  set  in  X-spacc. 


Lemma  1.1: 


/ 


A  point  is  an  extreme  point  of  the  polyhedral  set  in  X-space  if  and 


35 


only  if  it  corresponds  to  a  solution  with  a  single  flow-chain  between  each 
source-sink  pair. 


Proof : 


Suppose  there  are  n  source-sink  pairs  with  r^  >  0  .  Then  there 
are  n-equations  of  the  f.orn  (i.lb). 

Any  basic  feasible  solution  of  the  convex  polyhedron  defined  by 
equations  (1.3b)  and  (1.1c)  can  have  at  most  n  nonnegative  variables 
(because  tliere  are  n-equations  of  the  form  (1.1b)).  No  variable  appears 
in  more  than  one  equation  of  the  set  (1.1b)  and  all  the  n-equations  have 
to  be  satisfied.  So,  any  feasible  solution  must  have  at  least  n-positive 
variables.  Thus,  every  basic  feasible  solution  has  exactly  n-positivc 
variables,  one  for  each  source-sink  pair.  This  establishes  that  any  basic 
feasible  solution  corresponds  to  a  single  flow-chain  between  each  source- 
sink  pair. 

If  there  is  a  single  chain  between  each  source-sink  pair,  then  only 
one  variable  for  each  equation  (1.1b)  is  positive  and  each  of  these 
positive  variables  are  different  from  the  other.  So,  they  form  an 
independent  feasible  set.  Thus,  they  arc  a  basic  feasible  solution. 

So,  a  solution  is  a  basic  feasible  solution  if  and  only  if  it 
corresponds  to  the  solution  with  a  "'ingle  flow-chain  between  each  source- 
sink  pair.  A  solution  is  an  extreme  point  if  and  only  if  it  is  a  basic 
feasible  solution  (Page  100,  Reference  H-l).  Hence,  a  point  is  an 
extreme  point  of  the  polyhedral  set  if  and  only  if  it  corresponds  to  a 
solution  with  a  single  flow-chain  between  each  source-sink  pair. || 


In  this  problem,  a  concave  function  has  to  he  minimised  over  a 
convex  polyhedral  set.  The  following  lemma  substantially  reduces  the 
number  of  points  to  be  searched  in  any  direct  cimniorat ion  procedure. 
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Lemma  1.2: 


The  minimum  of  a  concave  function  on  the  convex  polyhedral  set  D  , 
if  it  exists,  is  attained  by  at  least  one  of  the  extreme  points  (vertices). 


Proof : 


Since  the  set  is  a  polyhedral  set,  any  point  Y  in  the  set  can  be 

expressed  as  a  convex  combination  of  its  extreme  points  Y, ,  ...»  Y 

l  m 

m 

Suppose  the  minimum  is  attained  at  Y  where  Y  =  £  a.Y  and 

o  o  i=1  l  i 

m 

]>  a.  =  1  ,  a  >  0  V  i.  Since  Y  is  the  minimum  point, 

i=l  1  1  ° 

f(Y  )  <  f (Y)  for  any  Y  c  D 
o  — 


(m  \  m 

I  a-iY-l  L  1  ajf(Y-)  because  f (Y)  is  a  concave 
i=l  1  7  i=l  1  1 

function. 

Let  f (Y  )  =  Min  f(Y.)  .  So, 
s  ±  i 


m  m 

f  (Y  )  >  £  a.f(Y.)  _>  l  a  f(Y  )  (since  a  0) 

°  i=i  1  1  i=l  1  s  x 


I  m 

*  f (Y  )  (since  ^  a  =  1 

8  \  i=l  1 


Since  Y  eD,f(Y)-f(Y).  Thus,  the  extreme  point  Y  must 
sos  s 

also  be  a  minimum  point. | | 


So,  a  direct  method  of  finding  the  minimum  of  the  concave  function  in 
a  polyhedral  set  is  to  enumerate  all  the  extreme  points  and  evaluate  the 
function  at  all  these  extreme  points  to  determine  the  minimum.  But  this 
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is  not  a  very  practical  method  because  the  number  of  extreme  points  in  the 
polyhedra  defined  by  (1.1b)  and  (1.1c)  is  very  large  and  grows  very 
rapidly  with  an  increase  in  the  number  of  arcs  and  source-sink  pairs. 

From  the  Lemma  1.2,  the  minimum  is  attained  at  one  of  the  extreme  points 
in  Y-space.  However,  the  extreme  points  in  Y-space  arc  very  difficult  to 
determine.  And  it  is  easy  to  determine  an  extreme  point  in  X-space  by 
using  a  single  chain  for  each  source  and  sink  pair  (from  Lemma  1.1). 
Therefore,  it  is  useful  to  determine  the  relationship  between  extreme 
points  in  Y-space  and  extreme  points  in  X-space,  which  is  done  in 
Lemma  1.3. 

The  image  of  a  point  P  of  X-space  in  Y-space  is  the  point  in  Y~ 
space  which  is  obtained  by  transforming  the  coordinates  of  the  point  P 
by  using  the  matrices  A_  . 

The  following  small  example  reveals  that  more  than  one  paint  in 
X-space  can  have  the  same  image  in  Y-space,  and  it  also  shows  that  the 
image  of  an  extreme  point  in  X-spacc  can  be  a  nonextreme  point  in 
Y-space. 


X 


T 


The  required  flows  are  r^ 

/  1357  13457  23456  2  )5b\ 

Vi 7  ,X17  ,X2G  ’*26  / 


*=  2  and  r«,  «  2  .  Let 
to 

and  vT  -  (yn.y23,y34,y35,y/t5,y56,y57) 


HUH 


18 


T 

The  vector  Y^  ,  where  =  (2, 2, 2, 2, 2, 2, 2)  ,  is  the  image  of  both 
and  ,  where  =  (2, 0,2,0)  and  X^  =  (0,2, 0,2)  .  Both  in  X^  and 
X^  there  is  a  single  path  between  each  source  and  sink  pair,  so  they  are 
extreme  points-  However,  Y^  is  a  convex  combination  of  two  feasible 
points  =  (2, 2, 4, 0,4, 2, 2)  and  Y^2  =  (2, 2, 0,4, 0,2, 2)  ,  since 

Y1  =  y  ^.1  +  \  1^2  *  So,  Y^  is  not  an  extreme  point. 

Lemma  1,3: 

For  a  point  to  be  an  extreme  point  of  the  polyhedron  in  the  Y-space, 
it  is  necessary,  but  not  sufficient,  that  it  be  the  image  of  an  extreme 
point  of  the  polyhedron  in  X-space. 

Proof : 

Suppose  the  lemma  is  not  true,  and  that-  there  exists  Y  which  is  an 

extreme  point  of  the  polyhedron  in  Y-space  and  it  is  an  image  of  X  e 

polyhedron  in  X-space,  which  is  not  an  extreme  point  (i.e.,  Y  is  the 

image  of  nonextreme  point  in  X-space).  So,  X  <=  £  a^X*  where 

i 

£  a  =  1  ,  a,  >.  0  V  i  ,  and  at  least  two  a.'s  are  nonzero,  and  X*’s 
i  1  1  1 

are  extreme  points  of  the  polyhedron  in  X-space.  Let  the  image  of  X*  in 

Y-space  by  Y^  .  Since  X  to  Y  is  linear  translation,  Y  =  £  a  Y*  • 

i  i 

Now  if  all  Y^’s  are  not  equal,  then  by  definition  Y  is  not  an  extreme 
point,  a  contradiction.  However,  if  they  are  all  the  same,  then  Y  is 
also  image  of  X*  which  is  an  extreme  point  in  X-space.  So,  the  necessary 
part  is  valid.  The  above  example  shows  that  it  is  not  sufficient. | | 

The  above  discussion  shows  that  if  we  se.trch  all  the  extreme  points 
in  X- space  wc  cover  all  the  extreme  points  in  Y-space.  Since  some  extreme 
points  in  X-spncc  may  correspond  to  nonext reme  points  in  Y-space,  this 
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search  procedure  may  involve  some  unnecessary  work,  but  the  ease  of 
generating  extreme  points  in  chain  space  justifies  the  emphasis  on  X- 
space  in  the  subsequent  algorithms. 

Another  major  difficulty  arises  in  pursuing  a  solution  to  the  above 
problem  due  to  the  fact  that  the  problem  has  quite  a  few  local  minimum 
points  which  are  not  necessarily  global  minimum  points. 

Definition: 

By  a  local  minimum  of  Z(Y)  on  D  (where  Y  e  Rn)  ,  we  shall  mean 
a  vector  Y  c  D  for  which  there  exists  a  positive  scalar  e  and  a 
corresponding  €-neighborhood  N(Y,e)  =  (Y  e  Rn  :  | | Y  —  Y | j  <  € }  ,  such 

that  Z (Y)  <  Z(Y)  for  all  Y  e  D  O  N(Y,e)  .  If  the  latter  inequality 
holds  for  all  Y  e  D  ,  then  Y  is  called  the  global  minimum  of  Z  in  D  . 

So,  any  gradient  search  method  which  looks  for  an  improvement  of  the 
objective  function  in  an  € -neighborhood  might  result  in  j  local  minimum 
point  (one  such  method  is  discussed  in  Chapter  3).  The  most  difficult 
aspect  of  this  is  the  identification  of  the  globally  minimum  point  when  it 
is  obtained.  Considerable  search  in  this  direction  convinces  the  author 
that  there  is  no  general  method  available  which  identifies  a  global 
optimal  point  for  such  optimization  problems,  though  such  identification 
in  some  special  cases  may  be  possible,  as  will  be  discussed  in  Chapter  2. 

Though  there  is  only  one  commodity,  telephone  calls,  the  problem  must 
be  treated  as  a  multi-commodity  problem  because  calls  are  distinguished  by 
their  starting  and  ending  points,  and  we  assume  t ha L  then-  art  many 
different  pairs  of  starting  and  ending  points  specified.  So,  the 
difficulties  of  multi-commodity  flow  problems  also  arise  here. 

If  (  as  in  the  synthesis  problem  we  are  studying)  there  is  no 
capacity  restriction  and  if  (as  is  not  the  case  in  our  problem)  the  cost 
is  a  linear  function  of  the  flow,  then  the  solution  can  he  easily  obtained 
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by  Floyd's  shortest  path  procedure.  Here  the  linear  cost  coefficients  are 
used  as  the  lengths  of  the  arcs  and  the  entire  flow  between  each  source- 
sink  pair  is  passed  through  the  shortest  path  between  them.  In  most  linear 
cost  problems,  however,  there  is  a  capacity  constraint.  In  such  cases,  a 
very  efficient  algorithm  exists  only  for  the  case  of  2  source-sink 
pairs  [H-3].  For  more  than  two  pairs,  a  certain  large  scale  programming 
method  using  a  column  generation  technique  is  useful.  This  linear  cost 
method  can  be  extended  successfully  to  convex  cost  functions  because  each 
arc  can  be  split  into  a  number  of  arcs,  which  transforms  the  problem  into 
a  linear  cost  problem  with  an  enlarged  number  of  arcs.  For  instance,  the 
convex  cost  on  an  arc  in  the  Figure  1.8  can  be  approximated  by  3  arcs 
having  linear  cost  coefficient  a^  ,  a^  ,  a^  ,  and  capacity  c^  *  °2  *  C3  ’ 
respectively. 


FTGURK  1.8 
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Here  the  proper  order  of  choice  of  arcs  is  guaranteed  (i.e., 
arc  (a1>Cl)  (cost  coefficient  and  capacity  c^  is  completely  filled 

before  (a2,c2>  is  used  and  when  (a2,c2>  is  filled  (a.^)  is  used) 
because  the  cost  coefficient  increases  with  flow  value.  However,  for 
concave  cost,  the  cost  coefficient  decreases  with  the  flow  value  so  the 
proper  order  will  not  be  maintained.  Hence,  such  simplification  is  not 
possible  and  no  procedure  of  finding  the  global  optimal  point,  which  does 
not  involve  the  risk  of  total  enumeration,  is  available.  The  methods 
proposed  in  succeeding  chapters  are  local  search,  and  branch  and  bound 
(here,  in  the  worst  possible  case,  total  enumeration  may  take  place) 
procedures. 
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CHAPTER  2 

The  last  chapter  exposed  the  practical  and  theoretical  difficulties 
involved  in  seeking  a  solution  procedure  for  the  problem  considered  in 
this  thesis.  Before  discussing  the  specific  solution  procedures  in 
succeeding  chapters,  we  will  survey  in  this  chapter  the  available  results. 
The  relevant  areas  of  research  which  must  be  covered  are  multi-commodity 
flow  problems  and  concave  cost  minimization  problems.  We  will  survey 
the  important  results  in  these  two  fields,  and  we  will  discuss  a  simple 
variation  of  this  problem  which  is  relatively  easy  to  solve. 

The  pioneering  work  in  multi-commodity  flow  problems  was  done  by,  among 
others,  L.  R.  Ford,  D.  R.  Fulkerson,  T.  C.  Hu,  R.  E.  Gomory  and 
J.  A.  Tomlin^. 

Ford  and  Fulkerson  [F-4]  initiated  the  analysis  of  the  maximum 
flow  problem  where  there  is  more  than  one  source  and  sink  pair.  They  showed 
that  if  the  flows  are  from  a  set  of  sources  to  a  set  of  sinks,  then  the 
maximum  sum  of  flows  between  the  two  sets  can  be  obtained  by  solving  a 
maximum  flow  problem  between  a  super  source  and  a  super  sink  in  an 
augmented  network,  where  there  are  extra  arcs  of  infinite  capacity  from 
the  super  source  to  the  set  of  all  sources  and  from  the  set  of  all  sinks 
to  the  super  sink.  They  also  indicated  [1-3}  a  solution  procedure  for 
the  more  general  problem  discussed  below,  gave  the  formulation  shown  in 
equations  (2.2a,  b,  c)  and  indicated  the  solution  method  discussed  there. 

The  more  general  multi-terminal,  multi-commodity,  maximum  flow 

problem  for  a  capacitated  network  was  formulated  by  Hu  and  Gomory  [G-4] 

as  follows.  Let  there  be  sources  N  and  sinks  N  ,(s  =  1,  ...,  q, 

s  s 

s’  =  1',  ....  q’)  where  flow  s  is  from  to  Ng ,  .  Let  be  the 

^The  names  of  R.  T.  Chien  [ C—  2 ]  and  W.  Maceda  [M-l]  also  should  be  mentioned 
in  this  connection.  But  for  various  reasons  1  will  not  be  discussing  their 
work. 
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flow  s  in  the  arc  ij  ,  and  f(s,s')  be  the  value  of  the  flow  s  from 

N  to  N  ,  •  The  first  type  of  problem  is  to 
s  s 

(2.1a)  Maximize  f(s,s') 

s=l 

if  j  =  s 

if  j  +  s  ,  s’ 

if  j  =  s’ 

(2.1c)  J  |xS  |  <b  (for  all  i,j) 

8=1  J  3 

(2. Id)  xS  0  (for  all  s,i,j)  . 


(2.1b) 


Subject  to 


l  xS. 

i  ^ 


I  Xs, 

k  Jk 


Unlike  the  one  commodity  max-flow  problem,  the  constraint  matrix  here 
2 

is  not  unimodular,  so  an  integral  solution  is  not  guaranteed.  Hu  [H-3J 
solved  a  special  case  of  the  above  problem,  the  case  of  two  commodities 

in  an  undirected  network  where  capacities  b„  are  even  integers  by 

I 

using  his  max  bi-flows  min-cut  theorem.  His  solution  procedure  involves 

the  following  steps.  First  the  labelling  algorithm  is  used  to  find  max-flow 

of  the  1st  commodity.  Then  a  cycle  of  flow  of  the  1st  commodity  is 

determined  in  such  a  way  that  by  introducing  some  flow  in  the  cycle,  the 

flow  values  in  certain  arcs  are  decreased,  so  that  a  flow  augmenting  chain 

of  the  2nd  commodity  is  obtained  to  introduce  the  2nd  flow  at  the  maximal 

possible  level  without  changing  the  flow  value  of  the  1st  commodity. 

Rothschild  and  Whinston  [R— 4 ]  have  extended  this  two  commodities  flow  result 

3 

for  a  pseudo-symmetric  network.  The  linear  programming  solution  of  the 


A  matrix  is  unimodular  if  all  the  subdeterminants  of  the  matrix  have  values 
0  or  1.  It  has  been  shown  (p.  125,  t H—  3 ) )  that  with  unimodular  matrix  D  . 
and  integral  vector  b  ,  the  convex  polyhedron  DX  >  b  has  all  integral 
extreme  points,  and  also  if  all  the  extreme  points  are  integral  and  b  is 
integral  vector  then  D  is  unimodular. 


3 


Definition: 
/ a  nod  e  i 


A  network  is 
ia  even  if 


Psuedo- lymmetric  if  it  has  nil 
£  b^  is  an  even  integerj. 


even  nodes 


max-flow  problem  in  the  multi- commodity  case  is  also  studied  by  Gomory 
and  Hu  [G-3] .  Here  the  arc-chain  formulation  is  used.  Consider  an  M 


arc  network  with  capacities  b, ,  . . . ,  b, ,  .  A  chain  in  the  network  can  be 

x  M 

represented  by  an  M  vector  with  1  in  a  component  if  the  arc  is  used, 
and  0  if  the  arc  is  not  used  in  the  chain.  Let  us  define  an  arc  chain 
incidence  m -trix  A  =  [a..]  as  follows 


a .  . 
ij 


1  if  the  arc  i  is  in  the  chain  j 
0  otherwise  . 


If  x.  Ts  the  amount  of  flow  in  chain  i 
J  J 

problem  is  given  by 


the  multi- commodity  max-flow 


(2.2a) 

(2.2b) 

(2.2c) 


Maximize  J  x. 

J  J 


Subject  to  £  aijXj  —  bi  d  =  ^  >  • • • >  M) 


x .  >  0 
J  “ 


where  b^  is  the  capacity  of  arc  i 


The  matrix  ^as  a  vcry  large  number  or  columns  (defined  as  a  ) , 

one  column  for  each  possible  chain  for  each  commodity.  Any  basis  has  only 
M  columns  and  at  each  step  one  has  to  consider  one  additional  column  as 
a  candidate  for  the  basis.  So  at  any  tine  we  need  to  consider  only  a 
(M  +  1)  *  (M  +  1)  matrix.  Suppose  we  have  M  columns  to  start  the 
algorithm  and  get  the  dual  price  vector  r  =  (it  ,  ...,  *  )  where  each 
Tij  corresponds  to  a  specific  row.  Notice  that  c  =  1  for  all  j  here. 
The  relative  cost  vector  of  every  noubnsic  column  a.  is  given  by 
c  -  1  -  Tin.  .  If  all  c.  <  0  then  the  present  basis  is  opt  intal.  If 
some  c.  >  0  wc  nerd  to  bring  that  column  into  the  basis.  All  of  the 
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tt^'s  can  be  made  nonnegative  by  introducing  into  the  basis  columns 
corresponding  to  slack  variables  in  rows  with  negative  n^'s  •  Then  the 
problem  of  determining  which  nonbasic  column  to  bring  into  the  basis  is 
simply  solved  by  the  column  generation  scheme  as  follows.  Use  n  as  the 
length  of  arc  i  and  find  the  shortest  chain  (a^)  connecting  any 
source-sink  pair.  If  all  na ^  >  1  ,  stop  since  all  c  <  0  and  we  have 
the  optimal  solution.  If  not,  bring  column  a.  into  the  basis. 

Hu  and  Gomory  also  considered  the  feasibility  problem  where  certain 

flow  requirements  have  to  be  satisfied,  which  Is  very  similar  to  the  previous 

problem.  Sometimes  the  problem  of  feasibility  is  coupled  with  the 

necessity  to  minimize  the  cost  of  building  the  capacity.  If  is  the 

cost  of  building  a  unit  capacity  in  arc  i  ,  we  would  like  to  minimize 

£  c.b^  •  This  problem  is  discussed  in  more  detail  by  Tomlin.  The  most 
i 

general  case  of  the  above  problem,  which  llu  solved,  is  the  case  in  which 
requirements  vary  over  different  time  periods.  Here  a  set  of  requirements 
has  to  be  satisfied  for  each  of  T  periods.  The  objective  is  to  minimize 
the  cost  of  building  a  sufficient  capacity.  An  example  of  multi-commodity 
minimum  cost  flow  problem  with  time  varying  requirements  is  solved  in 
Hu ' s  book  [H-3] . 

Tomlin  lT-1]  has  formulated  explicitly  the  problem  of  finding  minimum 
cost  flow  in  a  capacitated  network  which  satisfies  certain  flow  requirements 
for  a  directed  network  using  both  node-arc  and  arc-chain  matrices,  and  for 
an  undirected  network  using  only  the  arc-chain  matrix.  In  both  formulations 
he  obtains  a  large  scale  linear  program  having  a  special  structure  for  the 
constraint  matrix  so  that  the  Dantzig-Wol fe  decomposition  principle  can  he 
applied.  After  an  application  of  the  decomposition  principle,  a  subproblem 
(finding  the  shortest  route  problem  between  two  terminals)  must  be  solved 
In  order  to  generate  columns.  Many  efficient  algorithms  can  be  applied  to 
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solve  these  subproblems  because,  with  proper  manipulations,  we  can 
restrict  these  problems  so  that  they  have  only  nonnegative  arc  cost 
(shortest  route  problems  with  all  nonnegative  entries  in  distance  i-atrix 
are  easy  to  solve) . 

More  recently,  Richard  C.  Grinold  [G-5]  has  indicated  a  polyhedral 
game  approach  (similar  to  Kornai  and  Liptak  [K-3])  to  the  multi-commodity 
max-flow  problem  in  a  directed  network.  The  solution  procedure  involves 
allocating  flow  capacity  to  the  various  commodities,  solving  a  one- 
commodity  max-flow  problem  for  each  commodity  and  a  very  trivial  linear 
program  at  each  step.  The  method  is  easy  to  code  and  involves  simple 
computation  but  is  recommended  only  for  suboptimization  because  of  its 
poor  convergence  property.  He  also  indicated  an  extension  of  this 
procedure  for  the  multi-commodity  min-cost  flow  problem,  but  with  the 
restriction  that  in  order  to  make  the  problem  feasible  for  certain 
allocations  of  flow  capacity  we  need  to  introduce  high  cost  by-pass  arcs, 
parallel  to  the  original  arcs,  which  can  effectively  double  the  problem 
si  ze. 

So  most  of  the  work  in  multi-commodity  flow  is  restricted  to  max-flow 
or  to  linear  cost  cases  involving  a  capacitated  network.  The  problem 
considered  in  this  thesis  has  no  capacity  restrictions  -  so  a  major  set 
of  constraints  is  avoided,  thus  making  the  problem  comparatively  easier 
to  tackle.  However,  the  nonlinearity  of  cost  structure,  particularly 
concave  cost,  introduces  substantial  difficulties.  A  look  at  the  literature 
on  the  concave  cost  minimization  problem  is  appropriate  at  this  point. 

Philip  B.  Zwart  [ Z — 3 J  has  made  an  interesting  observation  for  a 
class  of  nonconvex  programming  problems  of  the  form: 

Minimize  F(X) 

Subject  to  G^(>.)  0 
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where  F(X),G  (X)'s  are  not  necessarily  convex  functions.  If  X  is  a 

local  minimum  and  A^'s  are  the  corresponding  lagrange  multipliers  for 

the  problem,  then,  if  the  modified  objective  function  F(X)  +  £ 

* 

becomes  a  convex  function,  X  is  also  a  global  minimum  for  the  problem. 

By  this  one  can  sometimes  recognize  when  a  local  minimum  point  is  a 
global  minimum  point.  Unfortunately,  if  we  do  not  have  convexity  then 
we  cannot  tell  whether  the  point  is  local  or  global  minimum. 

Willard  I.  Zangwill  [Z-l]  has  considered  a  special  type  of  flow 
problem  in  which  the  cost  of  flow  in  any  arc  is  a  concave  function  of 

■•V. 

the  amount  of  flow.  This  special  type  is  when  there  is  a  single  source 

and  a  number  of  sinks  or  a  single  sink  and  a  number  of  sources.  He  defined 

the  concept  of  extreme  flow  as  one  which  is  not  a  convex  combination  of 

two  other  flows.  The  extreme  flow  corresponds  to  the  extreme  point  of  the 

convex  polyhedron  generated  by  flow  conservation  conditions  in  node-arc 

or  arc-chain  space.  He  characterized  the  extreme  points  in  such  problems 

and  showed  that  extreme  flow  corresponds  to  arborescence  flow  in  multi-sink 

single  source  or  multi-source  single  sink  problems.  Based  on  this 

observation  he  developed  a  dynamic  programming  algorithm  to  solve  such 

problems.  Consider  a  network  with  source  node  1  and  sinks  a  and  b  , 

where  b  =  a  +  1  .  A  flow  from  node  1  of  r  units  to  a  and  r,  units 

a  b 

to  b  is  required  at  minimum  cost.  We  know  in  arborescence  flow  that 

some  arcs  will  have  r  +  r,  flow,  some  will  have  r  flow,  and  some 

a  b  ’  a 

will  have  r.  flow.  Furthermore  we  know  that  if  the  flow  r  +  r,  is 
b  a  b 

separated  at  a  node  into  rg  on  one  arc  and  r^  on  another,  the  flow 

will  never  be  r  +  r,  on  any  subsequent  arc.  Define  V  (a)  to  be  the 
ad  e 

minimum  cost  of  shipping  r  units  from  node  e  to  node  a  ,  V  (b)  to 

3  6 

be  the  minimum  cost  of  shipping  r^  units  from  lode  e  to  node  b  ,  and 


28 


V  (a,b)  be  the  minimum  cost  of  shipping  r  units  from  e  to  a  and 
e  a 

units  from  e  to  b  .  To  insure  that  no  material  is  shipped  from 
b  to  node  a  ,  define  V.  (a)  and  V,  (a,b)  as  very  large  numbers.  A(e) 

D  D 

is  defined  to  be  all  nodes  reachable  from  node  e  using  the  existing  arcs 
of  the  network.  (x)  represents  the  cost  of  x  amount  of  flow  in 

arc  (i,j)  ,  a  specified  function.  Then  the  recursive  relations  are: 


V  (a)  =  Min  {C  ,(r  )  +  V,(a)}  for  all  e  ^  a  ,  b  , 
e  ,  .  ,  .  efa  f 

feA(e) 


V  (b) 
e 


Min  {C  ,(r,)  +  V,(b)}  for  all  e  ^  b 
fcA(e)  ef  b  f 


V  <a,b)  =  Min  {C  f(r  +  r.)  +  V,(a,b), 
e  f,geA(e)  et  3  b  1 

C  r(r  )  +  V..(a)  +  C  (r,  )  +  V  (b) }  for  all  e  ^  a  ,  b 
ef  a  f  eg  b  g 


where  V  (a)  «  V.  (b)  -  0  ,  and  V  (a,b)  =  V  (b)  . 
a  b  a  a 

These  relations  are  solved  recursively  until  V^(a,b)  is  obtained, 

which  gives  the  minimum  cost  for  the  required  flow.  This  method  can  be 

generalized  for  more  than  2  destinations.  However,  with  n  sinks  the 

number  of  recursive  relations  to  be  evaluated  at  each  node  is  2°  -  1  . 

Hence,  the  calculations  required  become  prohibitive  for  large  n  . 

Arthur  F.  Veinott  [V-1,2]  has  considered  the  characterization  of  the 

4 

extreme  points  of  Leontief  substitution  system.  The  special  network  model 
considered  by  Zangwill  corresponds  to  a  transhipment  Leontief  substitution 
mode’..  So,  Veinott's  approach  puts  Zangwill's  algorithm  in  a  mere  general 
setting.  However,  except  for  special  cases  (e.g.,  the  one  considered  by 
Zangwill),  the  amount  of  computational  effort  required  to  search  the  extreme 
points  to  find  one  that  is  optimal  increases  exponentially  with  the  size 
of  the  problem  and  so  tends  to  be  enormous. 

4 

The  definition  of  Leontief  matrix,  as  well  as  transhipment  Leontief 
substitution  model  is  discussed  in  the  paper  by  A.  F.  Veinott  [V-1,2]. 


B.  Rothfarb  and  M.  Goldstein  [G-l]  in  a  recent  paper  considered 

the  one  terminal  Telpack  problem.  This  is  similar  to  Zangwill's  one 

sink,  multiple  source  network  problem,  but  they  considered  the  cost 

function  obtained  from  that  in  Figure  1.4  by  connecting  the  discrete 

points  by  straight  lines.  Assume  that  the  nodes  are  sequentially 

numbered  from  1  to  n  ,  where  the  node  n  is  the  sink  and  the  others  are 

the  sources.  If  an  arc  connects  nodes  i  and  j  ,  j  >  i  ,  it  will  be 

denoted  by  b(i,j)  .  If  the  flow  is  directed  toward  j  it  is  positive, 

and  if  directed  towards  i  it  will  be  negative.  The  cost  curves  are 

defined  by  means  of  their  points  of  discontinuity  in  derivative  (break 

points)  for  arc  b(i,j)  ;  let  [Wk(i,j) ,C^(i,j) ]  represent  the  flow  and 

cost  coordinates  of  the  kth  smallest  positive  value  of  flow  at  which  the 

cost  has  discontinuous  derivative.  Then  W *  -W^(i,j)  and 

C_k(i,j)  =  -Ck(i,j)  .  Let  CQ(i,j)  =  Wo(i,j)  =  0  .  The  level  of  flow  x 

in  an  arc  can  be  expressed  as  a  linear  combination  of  the  flow  levels 

at  the  nearest  break  point  above  x  and  the  nearest  one  below  x  .  To 

k=+K 

accomplish  this,  let  j)}  be  the  set  numbers  called  flow 

indices  associated  with  arc  b(i,j)  such  that 

+K 

0  =  -  1  and  £  “  1  » 
k  “  k=-K 

where  K  is  a  large  enough  number  to  cover  the  entire  required  range 

of  flow  values.  If  X  (i,j)  >  0  ,  then  only  X  1(i,j)  or  X  (i,j)  » 

but  not  both,  can  be  positive  and  only  these  two  indices  can  be  nonzero. 

K 

Then  the  flow  on  b(i,j)  can  be  represented  as  £  At,(i*j)Wi,(i*j)  • 

k=-K  k  * 

A  solution  of  this  problem  is  given  by  a  linear  program. 
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K 

Minimize  £  £  C  <i,j)-X  <i,j) 

(i,j)  k— K 

K  K 

Subject  to  £  £  W  <i,j)*X,  (i,j)  -  £  l  W  (j,i)*X  (j,i)  -  r. 

j|i<j  k=-K  k  *  j|i>j  k-K 

for  i  =  1,2,  ...»  n-1  , 

where  is  the  required  flow  from  the  soi’-ce  node  to  terminal  node  n  . 

Furthermore,  there  is  an  additional  constraint;  the  set  of  flow  indices 
k~K 

|X^(i,j)|  for  each  arc  should  satisfy  the  constraints  specified  above. 

k=-K 

So  the  ordinary  linear  programming  procedure  for  generating  the  solution 
is  not  enough.  Rothfarb  and  Goldstein  show  that  any  basis  of  the  linear 
program  will  have  arcs  of  two  kinds,  and  •  T^  are  the  arcs 

where  one  flow  index  is  1  and  the  rest  0,  and  are  the  arcs  with  2 

nonzero  flow  indices.  By  a  series  of  theorems  they  established  an 
intricate  procedure  which  treats  arcs  of  type  and  differently 

and  determine  which  arc  to  bring  in  the  basis  at  each  step.  So  the  pivot 
computation  is  far  more  elaborate  than  usual  simplex  method.  Also  the 
procedure  can  converge  to  a  local,  but  not  global  optimal  point.  Consequently, 
it  does  not  appear  superior  to  the  one-terminal  versions  of  the  algorithms 
of  Chapter  3  and  Chapter  4  to  follow. 

At  this  point  one  of  the  simple  variations  of  the  problem  which  is 
easily  solved  should  be  mentioned.^  In  this  case  there  is  just  a  single 
source  and  a  single  sink  with  flow  requirement  f  between  them  and  a 
concave  cost  function.  Since  it  has  been  established  that  there  must  be 
an  extreme  flow  which  is  the  minimum,  there  will  be  a  single  path  between 
the  source  and  the  sink.  The  flow  in  each  of  the  arcs  in  this  path  will 
re  f  .  To  find  the  optimal  path,  use  the  cost  of  f  amount  of  flow  in 

■*The  procedure  is  cryptically  mentioned  by  Zangwill. 

I 
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each  arc  of  the  network  as  the  length  of  that  arc.  The  shortest  chain 
between  the  source  and  the  sink  of  thl^>  network  will  be  the  global  minimum 
cost  solution. 

Hoang  Tui  [T-2]  has  described  a  procedure  that  seeks  the  global  minimum 
of  a  concave  objective  function  f(X)  with  the  polyhedral  constraint  set 
D  c  Rn  .  The  procedure  starts  with  a  local  minimum  point  X°  which  is 
assumed  to  be  a  nondegenerate  (a  necessary  condition)  extreme  point  of  the 
polyhedron  D  .  The  procedure  is  illustrated  for  n  c  2  case  in  Figure  2.1, 
where  D  is  OACEFG.  For  notational  simplicity  X°  (origin  0  in 
Figure  2.1)  is  taken  to  be  the  origin.  Along  its  n  distinct  edges,  n 


points  y 


....  y1,n  (y'1’*1  and  y^’^  along  OA  and  OG)  are  chosen 


lk  k  1  k  k 

such  that  y  *  =  £  *0  ’  where  £  is  the  direction  vector  of  the  kth 

edge,  6^’k  _  ^ax  |0  |  p(£ke)  >  f  =  f(X°)  ,  and  F(X)  is  the  concave 

extension**  of  f(x)  in  Rn  .  y^’^'s  are  linearly  independent  vectors. 
Within  the  cone  X°  ,  y^*\  y^,n  ,  (the  cone  Oy^’^y^'^),  the  value 

of  f (X)  ,  is  greater  than  or  equal  to  .  The  auxiliary  problem  at  this 
point  is  to  find  the  most  distant  extreme  point  of  the  polyhedron  D  from 


the  hyperplane  passing  through  y 


1,1 


,  y1,n  (the  line  y*’*y*’^  in 


Figure  2.1)  on  the  opposite  side  of  the  origin  X  (the  right-hand  side 
111° 

of  y  ’  y  *  in  Figure  2.1).  This  is  equivalent  to  solving  the  linear 
program 


n 

Maximize  h(X)  *=  ][  X 

k-1  k 


Subject  to  X  t  D 


Concave  extension  of  the  function  f(x)  is  any  function  concave  on  the 
whole  space  Rn  and  coinciding  with  f(x)  on  D  . 
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where  A^'s  arc  elements  of  the  n  dimensional  vector  A  =  B  *X  and 
B  is  a  matrix  with  columns  y^"*\  y^,n  . 

If  no  solution  of  the  above  linear  program  exists,  then  X°  is  the 
global  minimum  because  the  cone  defined  by  X°y^’\  ....  y~’1  contains 
the  entire  polyhedron.  Otherwise,  evaluate  f (X)  at  the  new  extreme 
point  (point  E  furthest  away  from  y^’^y^’^  in  Figure  2.1). 

Determine  =  Min  ja^,f(X^)|  (a^  <  in  our  example)  and  find 

2  k  2  12  2 

y  *  k  =  1,2,  ...,  n  along  n  edges  (Oy  ’  and  Oy  ’  in  Figure  2.1) 

2  k  -2  1 

as  before,  such  that  F(y  *  )  >  ,  and  also  determine  X  =  O^X  such 

that  0£  -  Max  |o  >  0  J  F(GX^)  >  a^j  .  Auxiliary  problems  involving 

lk  -2 

smaller  cones  arc  generated  by  replacing  one  of  the  y  *  by  X  and 

1  k 

keeping  the  others  fixed.  (Do  this  for  y  ’  only  in  case  the  corresponding 

Ak  *  0  .)  (For  the  problem  of  Figure  2.1,  the  two  auxil  iary  problems  are 

2  1-2 

finding  the  furthest  points  of  D  on  the  right  side  of  y  ’  X  and 
2  2-2 

y  ’  X  .  No  solution  exist  for  the  1st  problem  because  the  entire 

2  1-2 

polyhedron  is  on  the  left  of  y  ’  X  .  For  the  2nd  auxi  liary  problem 
2 

F(X  )  is  the  solution.  Now  use  point  F  in  the  same  way  that  E  was 
used  previously.)  In  this  way,  more  of  the  polyhedron  is  covered  by  each 
succeeding  set  of  auxiliary  problems;  when  the  entire  polyhedron  is  covered, 
we  have  obtained  the  global  optimal  solution.  The  advantage  of  this 
procedure  is  that  all  the  auxiliary  problems  have  the  same  constraint  set 
D  .  However,  the  procedure  is  not  useful  in  the  problem  of  this  thesis 
because  dimension  n  of  the  polyhedral  set  D  is  very  large.  The  number 
of  auxiliary  problems  generated  when  k  of  the  auxiliary  problems  arc- 
solved  is  1 11  ,  an  exponent  ial  growth.  Whereas  in  the  branch  and  bound 
procedure  des-ribed  in  Chapter  4,  the  number  of  sub-prohl eins  becomes  2 


when  the  kth  subproblem  is  solved.  Also  h(X)  is  given  in  such  a  way  that 
vc  car. not  use  any  efficient  method  like  the  shortest  path  method  to  solve 


the  auxiliary  problem. 


34 


K.  Ritter  [R-3,  C-3]  has  described  a  similar  method  of  determining 
the  global  optimal  solution  for  the  nonconvex  quadratic  cost  minimization 
problem  with  polyhedral  constraint  set.  The  procedure  involves  generating 
cutting  planes  which  exclude  local  optimal  points  without  excluding  the 
global  optimal  point  and  then  doing  a  local  search  (i.e.,  a  search  for  a 


local  optimal  point  in  the  reduced  set)  on  the  reduced  polyhedron.  This 
is  equivalent  to  a  variation  of  Hoang  Tui's  method,^  in  which  new  local 
search  procedure  is  carried  out  on  the  section  of  the  polyhedron  cut  out 
by  the  hyperplanes  generated  by  y^'\  ...»  y^,n  on  the  opposite  side 
of  origin.  The  generation  of  a  meaningful  cutting  plane  in  Ritter's 
method  itself  involves  solving  a  quadratic  programming  problem.  For  our 
problem  (assuming  quadratic  concave  cost  function) ,  both  the  generation  of 
cutting  planes  and  a  local  search  on  the  reduced  polyhedron  are  huge 
quadratic  programming  problems.  So,  the  procedure  is  not  that  useful. 
However,  as  we  will  point  out  in  Chapter  5,  some  clever  way  of  generating 
the  cuts  by  solving  simpler  subproblems,  and  also  a  simpler  local  search 
procedure,  may  be  a  hopeful  theoretical  direction  to  pursue  in  order  to  find 
the  global  optimal  solution  for  the  problem  of  this  thesis. 

A. Victor  Cabot  and  Richard  L.  Francis  [C-l]  have  described  a  method 
of  solution  for  the  nonconvex  quadratic  minimization  problem  by  ranking 
the  extreme  points.  The  method  per  se  is  not  applicable  for  our  problem, 
but  the  idea  can  be  adapted  to  yield  a  solution  procedure.  Consider  the 
problem: 


(PI) 


Minimize  f(Y)  -  CTY  +  YTDY 

Subject  to  Y  e  S  ,  where  S  =  (Y  |  AY  =  b,Y  >_  0} 


This  is  mentioned  in  the  paper  described  in  the  prec«_<_uing  paragraph. 
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where  Y  is  an  n-vector  of  variables,  the  matrix  A  is  m  *  n  ,  the 
matrix  D  is  n  *  n  ,  and  vectors  C  and  b  are  n-  and  m-vectors 
respectively.  d_.  j  =  i,2,  n  are  columns  of  the  D  matrix.  Solve 

u_.  =  Min  (dJ.Y)  ,  subject  to  Y  e  S  ,  and  assume  that  each  of  these  n 
problems  has  a  finite  minimum.  Then  consider  the  problem 


Minimize 


(P2) 


l  (c,  +  u  )y 

j=l  3  J  J 

Subject  to  Y  e  S 


where  y . ' s 
3 

vector  C  . 


are  elements  of  the  vector  Y  and 

^  /  T  \ 

Since  f(Y)  =  J  |c.y.  +  Y  d^y.j  , 


c.'s  are  elements  of 
3 

T 

and  Yd.  >  u.  ,  then 


g (Y)  <  f  (Y)  . 

If  Y°  is  the  solution  of  P2 ,  then  a  lower  bound  on  the  optimal 

value  f  of  PI  is  f^  =  g(Y°)  .and  an  upper  bound  is  =  f(Y°)  .  It 

k 

can  be  easily  proved  that  if  {Y  }  is  the  set  of  extreme  points  of  S 

k  * 

such  that  g(Y  )  <  f  then  the  optimum  solution  Y  of  PI  is  such  that 
=  u 

*  k-, 

Y  e  {Y  }  .  So,  the  algorithm  here  involves  searching  systematically  the 
set  (Y  }  .  To  find  the  extreme  points  of  S  yielding  the  2nd,  3rd, 

4th,  etc.  smallest  values  of  g(Y)  ,  Murty's  [M-4]  extreme  point  ranking 
method  is  used.  In  this  method  only  one  pivoting  step  (for  the  non- 
degenerate  case)  is  necessary  to  get  each  next  lower  ranked  extreme  point. 
Thus,  this  whole  procedure  involves: 


Step  0: 

Solve  n  linear  programs  to  obtain  u.  j  =  1,  ...,  n  . 


Step  1: 

Determine  Y°  ,  the  optimal  solution  of  P2.  Take  f  =  g(Y°)  and 
=  f(Y°)  and  Y°  as  the  current  best  solution  of  PI. 
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Step  2 : 

Use  the  extreme  point  ranking  procedure  to  get  the  next  best  extreme 
k  k 

point  solution  Y  to  P2.  If  g(Y  )  >  f  ,  then  stop  -  the  current  best 

*  k 

solution  is  the  minimum  solution  to  PI,  and  f  =f  .  If  g(Y  )  <  f  , 

u  ■=  u 

k 

then  replace  f  by  g(Y  )  . 

Step  3: 

k  k 

If  f(Y  )  <  f  ,  then  replace  f  by  f(Y  )  and  replace  the  current 

k 

best  solution  of  PI  by  Y  .  Otherwise  return  to  Step  2  without  changing 

f  or  the  current  best  solution, 
u 


Howev  ,  this  method  is  not  directly  applicable  to  the  problem  of  this 
thesis  for  the  following  reason.  The  objective  function  in  our  problem, 
when  a  quadratic  approximation  is  made  for  all  concave  functions,  is 


(Y)  -  £  (c.y.  +  d^2) 


where  d .  <  0 
J  = 


nd  y^  is  the  element  of  the  vector  Y  representing  flow  in  arc  j  . 
re  dj  is  a  scalar,  or,  in  the  notation  of  PI ,  D  is  a  diagonal  matrix, 
•nee,  u^.  =  Min  (d^y^)  subject  to  the  flow  restriction  of  our  problem, 
is  subproblem  is  to  find  the  minimum  cost  multicommodity  flow  in  the 
•twork  where  all  arc  costs  are  zero  except  one,  arc  j  ,  which  has  a 
incar  cost  with  negative  coefficient  d^  .  This  yields  a  negative  loop 
•o  the  objective  function  u^  tends  to  (by  passing  a  very  large  flow 

in  the  negative  loop) .  Thus  the  method  is  not  applicable  for  our  problem 
because  finite  u^  values  (which  are  unbounded  here)  are  necessary  to 
form  the  subproblem  P2.  However,  the  following  variation  of  the  procedure 
is  applicable. 


] 
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2.1  Solution  Procedure  by  Ranking  the  Extreme  Points 


The  problem  considered  has  the  objective  function 


M 

Minimize  F(Y)  =  £  f  (y . ) 

j=l  3  3 

Subject  to  Y  e  S 


where  S  is  the  feasible  polyhedron  in  arc  space  Y  i.e.,  Y  =  £ 

i.j 

?i>j 

A. .X..  and  X..  is  a  P..  vector  and  satisfies  the  constraint 

ij  ij  ij  l  ijf 

g 

equations  1.1b  and  1.1c.  If  we  can  find  the  upper  bound  Uj  of  the  flow 
in  each  arc,  then  we  consider  solving  the  problem 


M 

Minimize  G(Y)  =  C  y  where  C  U  *=  f(U.)  (in  Figure  2.2) 

j  2  J  J  J  J  J 

Subject  to  Y  e  S 


Then  G(Y)  <  F(Y)  in  the  entire  feasible  range  of  flow.  Solving  the 
modified  linear  problem  is  the  same  as  solving  shortest  chain  problems 
between  all  pairs  of  requirement  points  where  arc  lengths  are  C^'s  •  Let 
there  be  n  positive  requirement  (r^  >  pairs. 


Step  1: 

Find  the  shortest  chain,  using  C^'s  as  the  distances  of  the  arcs 

and  pass  the  required  flow  through  the  shortest  chain.  Let  Y°  be  the 

* 

arc  flow  vector.  Then  a  lower  bound  of  the  optimal  value  F  of  the  original 
problem  F  ’  C(Y°)  and  an  upper  bound  is  F  E  F(Y°)  .  The  current  best 
solution  is  Y°  . 


A  problem  faced  in  the  procedure  of  Chapter  4  and  discussed  there.  A 
trivial  upper  bound  is  the  sum  of  all  required  flows. 
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Determine  the  2nd  shortest  chains  between  all  n  pairs  of  requirement 

points.  Let  A^  ,  be  the  difference  between  the  2nd  shortest  chain  and 

the  shortest  chain  between  source-sink  pair  ij  .  The  list  of  n  numbers 

(A.  ,r. .}  is  the  amount  of  increment  for  the  neighboring  extreme  points  of 
ij  ij 

Y°  ,  one  (the  smallest)  of  which  corresponds  to  the  increment  of  the  next 
best  extreme  point  of  S  . 


Step  2: 

Let  s  and  t  be  the  source-sink  pair  which  gives  the  smallest 

increment  (from  the  list)  over  Y°  ,  among  all  neighboring  extreme  points 

o  k-1 

of  the  extreme  points  (Y  ,  . . . ,  Y  )  so  far  determined.  And  let  it  be 

the  neighbor  of  Ym(m  <  k  -  1)  .  Then  displace  the  flow  r  from  the 

chain  between  s  and  t  corresponding  to  the  point  Y™  to  the  newly 

found  chain  between  s  -  t  which  gives  minimum  increment.  Calculate  the 

k  k 

new  arc  flow  vector  Y  .  If  G(Y  )  >  stop  -  the  current  best  solution 

* 

is  the  optimal  solution  and  F  *  F  . 


Step  3: 

Solve  the  next  best  chain  between  s  and  t  ;  find  A  the 

st 

difference  between  previously  found  chain  and  currently  found  chain 

k 

between  s  and  t  .  One  neighboring  extreme  point  of  Y  has  increment 
(G(Yk)  -  G(Y°)  +  A  *r  \  over  Y°  and  remaining  (n  -  1)  neighboring 
extreme  points  have  increments  (G(Yk)  -  G(Ym)  +  increments  over  Y°  of 
neighboring  extreme  points  of  Ym  except  Yk}  over  Y°  .  Include  all 
these  increments  in  the  list. 

Step  4: 

If  C(Yk)  <  Fu  then  replace  F^  by  G(Yk)  .  If  F(Yk)  <  Fy  replace 
F^  by  F(Y  )  and  current  best  solution  of  the  original  problem  by  Y 
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Otherwise  go  to  Step  2  without  changing  or  the  current  best  solution. 

So,  in  this  procedure,  after  finding  the  2nd  best  chain  for  all 

pairs,  we  need  to  find  only  the  next  best  chain  for  a  particular  node 

pair  at  each  iteration.  To  reduce  the  number  of  searches  we  could 

simply  concentrate  only  on  loopless  chains.  A  chain  with  a  loop  cannot 

possibly  be  a  minimum  for  the  original  problem  because  by  taking  out  the 

flow  in  the  loop  we  can  reduce  the  cost.  J.  Y.  Yen  [Y-2]  has  proposed  a 

good  method  of  finding  the  K  best  loopless  chains.  However,  to  avoid 

confusion  one  should  note  that  a  chain  with  a  loop,  which  is  an  extreme 

point  in  arc  chain  space,  can  also  be  an  extreme  point  in  node  arc  space. 

So  by  ignoring  looped  chains,  we  may  be  ignoring  some  extreme  points, 

but  these  are  obviously  uninteresting  ones.  At  any  stage  k  in  this 

procedure  we  need  to  store  only  increments  of  the  neighloring  extreme 

points  of  Y^  over  Y°  such  that  G(Y°)  +  (increment  of  neighboring 

Ic  o 

extreme  points  of  Y  over  Y  )  <  F  and  no  others.  The  relative  merit 

c  U 

of  this  method  compared  to  the  methods  described  in  succeeding  chapters 
depends  on  the  efficiency  of  the  kth  shortest  chain  determination 
procedure.  For  the  general  problem  such  comparisons  are  hard  to  make. 

As  a  last  note  on  the  available  literature  on  this  TELPAK-type  problem, 
we  should  mention  the  work  of  G.  C.  Watling  and  J.  H.  Weber  [W-2].  They 
developed  a  heuristic  algorithm  which  synthesizes,  from  the  traffic  data, 
the  network  and  also  the  best  position  for  switching  centers.  This 
heuristic  procedure  uses  only  the  total  amount  of  originating  and  terminating 
traffic  as  the  input  traffic  data  and  does  not  consider  the  traffic  flow 
between  specific  originating  and  terminating  points.  This  is  reasonable 
if  the  traffic  is  statistically  well  dispersed  between  all  points,  and 
this  method  can  handle  a  very  large  problem.  It  is  purely  heuristic  and  can 
yield  a  solution  that  is  not  even  an  extreme  point  of  the  problem  we  are 
studying,  and  is  not  very  mathematical  in  approach. 
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CHAPTER  3 

The  particular  solution  procedure^-  described  in  this  chapter  can  be 
used  to  solve  moderate  sized  problems  (networks  of  up  to  200  nodes)  and 
it  is  very  efficient  (i.e.,  the  procedure  generally  converges  in  a  few 
iterations).  However,  it  yields  locally  optimal  points  and  one  does  not 
know  how  near  the  value  is  to  the  global  optimal  value.  In  this  chapter 
strategies  that  are  useful  in  obtaining  the  global  optimal,  using  a 
procedure  that  yields  local  optimal  points,  are  described.  Certain 
statistical  procedures  will  also  be  developed,  including  a  sequential 
sampling  plan  where  further  search  for  the  global  optimal  point  is 
stopped  when  the  cost  of  further  computation  becomes  more  than  the 
estimated  gain  in  reduction  in  the  optimal  value. 

For  the  main  procedure,  of  this  chapter,  to  work,  all  the  functions 

f  (y  )  ,  the  cost  of  y  amount  of  flow  in  arc  m  ,  have  lo  be 
m  •'m  ’  •'m 

nondecreasing  concave  functions  in  the  feasible  regions  of  Ym  •  In 
Appendix  D  it  is  shown  that  any  concave  function  is  continuous  everywhere 
except  perhaps  at  the  left  and  right-hand  end  points  of  the  interval. 
Since  the  function  is  nondecreasing  it  has  to  be  continuous  at  the  right- 
hand  end  point.  This  continuity  property  may  not  hold  at  the  left-hand 
end  point;  there  could  be  a  jump  at  y^  =  0  .  Appendix  D  also  proves 
that  both  left-hand  and  right-hand  derivatives  D  f  (y  )  and  D+f  (y  ) 
exist  at  all  points  except  the  end  points  of  the  interval,  and  the 
following  relations  hold: 


This  procedure  was  suggested  by  Bernard  Ynged  [ Y— 1 ]  fur  concave  cost 
functions  which  have  both  first  and  second  derivatives,  hut  we  shall  show 
that  it  is  valid  for  any  nondecrcas ing  concave  cost  function. 
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(3.1a) 

(i) 

V 

2 

^m 

>  y1 
m 

(3.1b) 

(ii) 

°+ 

i-h 

a 

ej- 

II  V 

« 

H>+ 

V 

2 

1 

>  y 
'm 

(3.1c) 

(iii) 

D"f  (yj  >  D+f(y) 
m  m  =  mm 

V 

ym 

in  the  open  feasible 
interval  . 

The  nondecreasing  property  of  the  function  guarantees  that  both  left 

and  right-hand  derivatives  are  nonnegative.  The  following  discussions 

and  propositions  characterize  the  optimal  flow  pattern. 

By  Lemmas  1.1  and  1.2  of  Chapter  1  we  have  shown  that,  if  the  cost 

functions  are  concave,  in  the  optimal  flow  pattern  there  is  a  single 

chain  between  each  source-sink  pair;  i.e.,  for  every  source-sink  pair 

* 

st  ,  a  required  flow  r  passes  through  a  single  chain  p  in  the 

St  St 

minimum  cost  (optimal)  network.  (However,  degeneracy  may  occur;  two 
chains  may  have  the  same  cost.)  We  will  show  a  stronger  condition  for 
flow  patterns  where  cost  functions  are  concave. 

Proposition  3.1: 

If  there  are  two  chains  carrying  total  flow  r  between  a  single 

source-sink  pair  st  then  the  entire  flow  r  can  be  put  in  one  of 

st 

the  two  chains  without  increasing  the  total  cost. 

Proof: 

Suppose  there  are  two  chains  P  and  P,  between  source-sink  pair 

a  b 

st  carrying  and  x^  amount  of  flow  respectively  such  that 

+  ■  r.t  • 

For  any  arc  m  in  chain  P  with  flow  level  y  +  x  ,  the 

m  m  a 

contribution  to  the  cost  due  to  flow  between  s  and  t  in  this  arc  is 

taken  to  be  g_(x  )  ■  f  (y  +  x  )  -  f  (y  )  .  A  similar  definition  exists 
m  a  m  m  a  m  m 
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for  the  arcs  in  chain  .  g^C* )  is  also  a  concave  nondecreasing 

function  because  g  (•)  is  obtained  from  f  (•)  by  shifting  the  origin 

m  m 

to  y  .  Define  the  functions  C_.  (x  )  =  £  g_(x  )  and 

'm  P  a  a 

a  mtP 

a 

Cp  (x^)  *  ^  g^x^)  '  ^rom  above  definitions,  the  cost  of  r 

b  meP. 

D 

flow  is  taken  to  be  Cp  (xa)  +  Cp  (x^)  and  both  Cp  (•)  and  Cp  (•) 

a  b  a  rb 

are  positive  sums  of  concave  nondecreasing  functions,  hence  they  are  also 
concave  and  nondecreasing.  Therefore, 


X  X 

C  (x  )  >  — -C  (r  J  +  —C  (0) 
P  a  **  r  P  st  r  P 
a  st  a  st  a 


cp.  (V  i  rrS  <rst>  +  rrcp<0)  because  rst  -  +  s 

b  st  b  st  b 


(3.2) 


Cp  fra)  +  C  <xb)  >  — .cp  (r3C)  +  (rst) 

a  b  st  a  st  b 


since  C  (0)  and  C  (0)  are  zero.  The  scalars  C  (r  )  and  C  (r  ) 

f  _  Jr.  *  St  r.  St 

a  b  a  b 

can  have  one  of  the  following  relations: 


(i)  S  <rst>  •  S  <v> 

a  b 


(il)  Cp  (r  )  >  C  (r  > 
a  b 


(111)  Cp  (rst)  <  C  (rat> 
a  b 


If  (i)  occurs  then  putting  the  entire  r  flow  in  P  or  P.  results 

st  a  b 

in  no  cost  increase  by  the  inequality  (3.2) • 
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If  (ii)  occurs,  then  putting  the  entire  flow  in  P.  results  in  less  cost 

D 

because 


.  .  *b 


CP  '*.>  +  S  <*b>  i  T~Ct  <'.t>  +  r~'CP  <r.t>  >  CP  <rst> 
a  b  st  a  st  b  b 


And  if  (iii)  occurs,  then  putting  the  entire  flow  in  P&  results  in  less 
cost  because 


cp  (*a> +  W  i  rr'cP  (rSt> +  r;-cp  -« 

■  st 


(r8t>  *  CP  (rst> 
a  a 


Thus  the  entire  r  .  flow  car  I»e  put  into  one  of  the  chains  without 
8t 

any  cost  Increase. | | 

Corollary  3.1: 

If  there  are  multiple  (more  than  one)  chains  between  a  source-sink 
pair,  they  can  be  replaced  by  a  single  chain  (one  of  the  multiple  chains) 
without  increasing  the  total  cost. 


Notation : 


pflb  c  p^  means  the  arcs  of  the  chain  pab  between  a  and  b  is 
a  subset  of  the  arcs  of  the  chain  p^  between  i  and  j  . 


Proposition  3.2: 

In  the  optimal  flow  pattern,  if  the  optimal  chain  p*^  between 

source  1  and  sink  j  passes  through  nodes  a  and  b  ,  and  there  is  also 

a  flow  requirement  between  a  and  b  ,  then  the  optimal  chain  p* 

ab 

,  *  * 
between  a  and  b  is  such  that  p  .  C  p. ,  . 

rab  rtJ 
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Proof : 

Let  r  ^  be  the  required  flow  between  i  and  j  and  r^ 
be  the  required  flow  between  a  and  b  .  Suppose  a  precedes  b 
in  the  chain  from  i  to  j  .  The  flow  pattern  is  unchanged  if  we 


consider  requirements  r^a  *  r^  , 

r *  *  r  +  r  ,  ,  and 

ab  ij  ab 

r 1  *  r 

bj  rij 

• 

Since  there  is  a  single  chain  between  i 

and 

a  ,  a  and 

b  ,  and 

b 

and  j  ,  chain  ab  coincides  with 

chain 

ij 

between  a 

and  b  , 

and 

*  * 

thus  p^  c-  P^j  •  if  b  precedes 

a  in 

the 

chain  from 

i  to  j 

then 

flow  pattern  is  unchanged  if  flows 

r  ’  - 

ib 

"ij 

*  rab  *  rij 

+  r  ,  ,  and 

ab 

r'.  «*  r..  .  Since  there  is  a  single  chain  between  i  and  b  ,  b  and 
aj  ij 

a  ,  and  a  and  j  ,  Pflb  C  p  . | | 


Thus  we  can  restrict  our  search  f jr  the  optimal  flow  network  to 
solutions  which  have  the  following  two  properties:  (i)  a  single  chain 
between  each  source  and  sink  pair,  and  (ii)  p^  c  if  a  source  and 

sink  pair  ab  is  contained  in  the  chain  of  another  source  and  sink  pair 

ij  • 


3.1  Definition  of  g-Optimal  Routing 

A  routing  (flow  pattern)  is  called  €-optimnl  if  it  has  the  above 
two  properties,  and  if  the  least  cost  chain  for  an  additional  flow  of 
€  units  (where  €  is  an  arbitrarily  small  positive  real  number)  between 
source-sink  pair  ij  is  the  same  as  the  chain  taken  by  the  ij  flow  of 
r^j  units. 

The  following  example  will  clarify  the  concept  of  c -optimality. 

For  the  network  shown  in  Figure  3.1b,  the  costs  of  chains  A  and  B  are 
shown  in  Figure  3.1a  by  curves  A  and  B  respectively.  R  is  the  point 
where  the  tangent  of  curve  B  has  the  same  slope  as  A  .  Suppose 
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r  =  R,  <  R  and  the  entire  flow  is  in  chain  B  .  For  an  additional 
st  1  = 

flow  of  amount  C  the  cheapest  chain  is  A  because  the  slope  at  0 
of  A  is  smaller  than  the  slope  at  of  B  .  Thus  chain  B  is  not 

c-optimal.  But  if  the  entire  flow  is  in  A  then  an  additional  flow  of 
amount  e  will  take  chain  A  because  the  slope  at  0  of  B  is  greater 
than  the  slope  at  R^  of  A  ,  and  chain  A  is  e-optimal.  It  is  also 
globally  optimal.  Suppose  r  =  R£  >  R  ,  and  the  entire  flow  is  in  B  . 

An  additional  flow  of  amount  6  will  take  path  B  because  the  slope  at 
R^  of  B  is  smaller  than  the  slope  at  0  of  A  .  Similarly,  if  the 
entire  flow  is  in  A  ,  then  an  additional  flow  of  amount  e  will  take 
chain  A  because  the  slope  at  0  of  chain  B  is  greater  than  the 
slope  at  of  chain  A  .  So,  both  A  and  B  are  e-optimal  chains. 

However,  only  A  is  globally  optimal. 

The  algorithm  discussed  below  searches  for  an  e-optima]  point.  (The 
relationship  between  e-optimality  and  local  optimality  is  discussed  later.) 
The  concept  of  e-optimality  reveals  a  clue  to  the  constructive  approach 
(of  the  algorithm)  which  the  conventional  definition  of  the  local  optimum 
does  not  provide. 

The  following  notations  are  used  in  the  subsequent  discussions.  P 

is  an  M- vector  where  the  element  p  is  the  amount  of  flow  in  arc  m 

m 

corresponding  to  a  feasible  flow  pattern  (i.e.,  a  pattern  satisfying  (1.1) 
of  Chapter  1).  e  is  an  M-vcctor  of  0's  and  l's.  If  e  -  1  indicates 
the  inclusion  of  arc  m  and  em=  0  indicates  its  exclusion,  then  any 
chain  between  any  source-sink  pair  can  be  represented  by  a  vector  e  . 

C(P)  is  the  total  cost  of  the  flow  on  all  arcs  corresponding  to  point  P 

/  M  \  . 

(i.e.,  C(P)  =  7  f  (p  )}  .  6  is  a  scalar  and  (P  +  6-e1  -  5*e  )  is  a 

\  m=l  / 

point  where  6  amount  of  flow  is  taken  out  from  chain  c 


of  point  P 


and  redirected  to  chain  e"^  .  F  (P)  represents  the  total  cost  of  all 

the  arcs  in  chain  e  of  point  P  ji.e.,  £  f  (p  ) ]  ,  and  it  is 

\  racchain  e  / 

a  concave  function  of  P  . 

The  following  small  example  of  Figure  3.2  shows  a  point  can  be 

C-optimal  but  not  local  optimal.  Figure  3.2a  gives  the  cost  function 

f .(•)  and  f  (•)  for  the  flow  in  chains  A  and  B  respectively  of  the 

network  3.2b  between  source-sink  pair  st  with  flow  requirements  r  ^  . 

Suppose  the  entire  flow  is  in  chain  B  .  And  the  left  and  right-hand 

derivatives  of  cost  functions  in  their  present  flow  levels  have  the 

following  relation:  D  £  (r  )  >  D+f  (0)  >  D+f  (r  )  .  Then  an  additional 

DSC  A  d  St 

flow  amount  e  will  take  the  chain  B  because  D+f.(0)  >  D+fn(r  )  . 

A  B  st 

Hence  the  present  flow  pattern  is  c-optimal.  But  if  c  amount  of  flow 

is  displaced  from  B  to  A  then  total  cost  decreases  because 

c*D  f  (r  )  >  e*D+f„(0)  .  So,  it  is  not  locaJly  optimal  point. 

B  st  A 

Lemma  3.1: 

For  a  concave  cost  function  a  local  optimal  point  having  a  single 

chain  between  each  source-sink  pair  is  an  e-optimal  point.  An  e-optimal 

* 

point  is  a  local  optimal  point  only  if  the  single  path  e  between  every 
source-sink  pair  also  satisfies  the  following  inequality 

(3.3)  D_F.(P)  <  D+F  (P)  V  e 

e*  =  e 

where  e  connects  the  same  source-sink  pair  as  e  .  (This  inequality 
is  satisfied  by  any  e-optimal  point  for  problems  with  concave  cost  function 
where  derivatives  exist  at  all  points.  Hence  in  this  case  a  point  is  a 
local  optimum  if  and  onl>  if  it  is  c-optimal.) 
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Proof : 

Consider  a  locally  optimal  point  P  having  a  single  chain  between 

each  source-sink  pair.  Consider  a  point  Q  within  e-neighborhood  of  P  . 

Q  is  obtained  from  P  by  redirecting  some  6(<  e)  amount  of  flow 

k 

between  a  source-sink  pair  from  the  single  chain  e  to  some  other  chain 

e  .  (There  could  be  more  than  one  chain  and  could  be  more  than  one 

source  and  sink  pairs.  But  these  cases  of  more  than  one  different  chain 

can  be  reduced  to  only  one  different  chain  between  a  source-sink  pair 

k 

without  any  cost  increase  by  Corollary  3.1.)  Q  =  P  +  6-e  -  6-e  .  Let 

*  * 

=  e  -  e  He,  and  F_*  and  F  A  be  defined  as  cost  of  flows  in 

c  e  He 

—  k  k 

the  set  of  arcs  in  e  and  e  fl  e  respectively. 


C(P  +  6-e  )  =  C(P)  +  F  .(P  +  6-e  )  -  F  .(P) 


C(P)  +  F_*(P  +  5-e*)  +  F  *  (P  +  6-e*)  -  F_*(P) 
e  e  He  e 


F  *  (P) 

e  ile 


C(P  +  6-e)  = 


k  k  k 

C (P  +  6-e  -  6-e  +  6-e  )  =  C(Q  +  6-e  ) 


=  C (Q)  +  F_*(Q  +  6-e*)  +  F  *  (Q  +  6-e*)  -  F_*(Q) 
e  e  fie  e 


F  *  (Q) 


=  C(Q)  +  F_*(P)  +  F  *  (P  +  6-e*)  -  F_* (P  -  6-e*) 
e  e  fie  e 


F  *  (P) 

e  fie 


The  cost  value  F_^.(Q  -1-  6-e  )  depends  on  cost  of  flow  in  the  arcs  in  the 
e 

—  k 

set  e  .  At  Q  the  flows  are  reduced  from  the  level  at  P  by  6  on 
—  * 

this  set  e  ,  and  when  6  is  added  back  the  level  of  flow  in  these  arcs 

* 

is  the  same  as  at  P  ,  so  F  ^(Q  +  6-e  )  =  F  ^ (P)  .  Similarly,  we  can 

e  e 

show  F  ^  (Q+6-e)=F<c  (P+6-e)  ,F  *(Q)  =FJk(P-6-e)  ,  and 

e  lie  e  He  e  e 

F  *  (Q)  =  F  *  (?)  . 

e  fie  e  fie 
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Now,  C(P  +  6* e)  -  C(P  +  6-e*)  =  C(Q)  -  C(P)  -  2/F  *(P)  - 

'  c 

%’ F_^ (P  +  6-e  )  -  h' F* (P  -  6-c  )\  >  C(Q)  -  C(P)  (because  F_A(-)  is  a 
ee  e 

concave  function)  >  0  ,  because  P  is  a  local  optimal  point.  So, 


C(P  +  6-e  )  <  C(P  +  6-e) 


or. 


{C(P  +  6-e*)  -  C(P)}  <  {C(P  +  6-e)  -  C (P) }  . 

Hence,  extra  6(<  €)  flow  is  passed  through  the  optimal  single  chain. 

So,  P  is  an  €-optimal  point. 

Let  the  e-optimal  point  P  satisfy  the  inequality  (3.3).  Consider 
a  point  Q  in  e -neighborhood  of  P  such  that  Q=P+6-c~6-e  where 
6(<  O  . 

C(Q)  =  C(P  +  6-e  -  6 • e*)  =  C(P)  +  6-D+F  (P)  -  6-D~F  /P)  >  C(P) 

e  e?‘  — 

(because  (3.3)  is  satisfied).  Any  point  Q  in  this  neighborhood  will 
have  the  above  relationship,  so  P  is  a  local  optimal  point. 

A  point  P  is  €-optimal  if  the  following  inequality  is  true  for  V  e 

(C(P  +  6 •  c* )  -  C(P)}  <  {C(P  +  6-e)  -  C(P>> 


or, 

6-D+F  ^ (P)  <  6-D+F  (P)  /because  C(P  +  6-e*)  =  C(P)  +  6-D+F  *(P) 
e  \  e 

or,  D+F  A(P)  <  D+F  (P)  V  e  .  If  the  derivative  exists  at  all  points 
e 

-  +  + 

then  D  F  ^(P)  =  D  F  ^ (P)  <  D  F  (P)  .  So,  the  6-optimal  point  also 
e  e 
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satisfies  the  inequality  (3.3),  and  therefore  is  a  local  optimal  point. 
Thus  for  this  case  a  point  is  e-optimal  if  and  only  if  it  is  a  local 
optimal  point. 

The  example  of  Figure  3.2  shows  a  counterexample  if  the  inequality 
does  not  hold. | | 

3.2  Algorithm 

In  the  description  of  the  algorithm  the  derivative  Df(«)  of  the 
cost  function  at  a  point  is  used.  If  this  derivative  does  not  exist  then 
the  right-hand  derivative  D+f(*)  is  used. 

Step  0; 

An  arbitrary  positive  number  representing  length  is  assigned  to 
each  arc  of  the  network  (e.g..  Euclidean  distance.  Or  an  arbitrary 
flow  level  is  assigned  for  each  arc  and  the  derivative  of  the  cost 
function  at  the  flow  level  is  used  as  the  length  of  the  arc). 

Step  1: 

Using  the  lengths  specified,  the  shortest  chain  is  determined 
between  each  source-sink  pair  (if  there  are  many  pairs  of  sources  and 
sinks,  Floyd's  algorithm  is  used). 

Step  2: 

The  entire  required  flow  between  a  source-sink  pair  is  passed 
through  the  shortest  chain  between  them.  The  total  flow  y  in  each 
arc  m  is  determined.  The  total  cost  of  the  flow  is  determined.  If  the 
total  cost  is  unchanged  from  the  previous  iteration  t lie  procedure  is 
stopped.  (Note:  If  the  flow  in  each  arc  remains  unchanged  in  two 
consecutive  iterations,  then  the  shortest  chain  and  the  total  cost  also 
remain  the  same.  However,  the  flow  might  change  between  two  degenerate 
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points  without  having  any  change  in  the.  total  cost.  Cycling  may  occur 
if  we  use  only  the  flow  value  to  find  the  termination  point.  Such  cycling 
is  not  possible  if  the  cost  function  is  strictly  concave.) 


Step  3: 


Derivatives  of  the  cost  function,  Df  (y  )  ,  at  the  flow  values  y  , 

m  -'m 

are  determined  for  each  arc.  D+f  (y  )  for  functions  having  no  derivative 

m  ro 

at  y^  .  These  derivatives  are  used  as  the  new  lengths  of  the  arcs. 


Go  to  Step  1. 


3.3  Remarks 


When  the  algorithm  stops,  then  the  specified  chain  p„  between  a 

source-sink  pair  (ij)  satisfies  the  inequality  [  Df  (y  )  < 

^  m  m  = 

mep .  . 

£  Df  (y  )  V  p..  chain,  assuming  derivatives  exist.  Now, 

mep .  .  J 

ij 

=  D+fm(ym)  =  D  fm(y •  The  inequality  £  D+f  (y  )  < 
m  'm  mm  m  m  m  'm  = 

mEpij 

)  D  f  (y  )  is  true  even  if  Df  (y  )  does  not  exist.  Therefore,  the 
^  m  m  m  m 

mep 

3  /  + 
point  is  an  ^-optimal  point  J f rom  the  proof  of  Lemma  3.1,  where  D  F  *(P) 

\  e 

is  the  same  as  £  D+f  (y  )  and  y  represents  point  P  and 

^  m  m  m 

mePij 

*  *  \ 

p. .  =  e  I  .  The  inequality  (3.1)  also  is  satisfied  if  Df  (y  )  exists: 

ij  / 

in  that  case  the  point  is  a  local  optimum.  However,  a  post-optimization 

procedure  is  necessary  to  obtain  a  local  optimal  point  for  functions  whose 

derivatives  do  not  exist.  Take  each  source-sink  pair  and  change  the 

lengths  of  arcs  along  the  optimal  single  chain  from  D+f  (y  )  to  D  f  (y  ) 

mm  mm 

If  the  shortest  route  remains  unchanged,  then  (3.3)  is  satisfied  and  the 


point  is  a  local  optimum.  If  the  shortest  chain  changes,  then  an  improvement 
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(reduction)  in  total  cost  can  be  achieved  by  using  this  new  shortest  chain. 
The  steps  of  the  algorithm  are  repeated  using  this  new  point.  Cycling 
is  not  possible  here,  because,  if  a  different  shortest  chain  is  obtained  in 
this  procedure,  the  total  cost  strictly  decreases  when  the  flow  is  redirected. 

The  motivation  for  using  derivatives  of  the  cost  function  is  the 
following  Kuhn-Tucker  necessary  condition  for  optimality  for  differentiable 
cost  function  (i.e.,  Df(*)  exists).  The  unrestricted  minimization 

problem,  equivalent  to  the  problem  (1.1)  with  dual  variables  X  for 
constraints  (1.1b)  and  8  for  nonnegativity  constraints  (1.1c)  is 
given  by: 


L  ■ "  W 

m=l 


l 

All  ij 


ij 


k=l 


ij 


-  r 


ij, 


i  ¥ 

All  ij  k=l 
>i>j 


k  k 
3. ,x .  . 
ij  iJ 


Since  each  variable  y 

m 

coefficients  0  or  1, 


is  a  linear  combination  of  x„ 


variables  wi th 


3L 


9x. 


I  DW 

ntpij 


-  X 


ij 


3k.  =  0 
ij 


The  above  equation  and  the  Kuhn-Tucker  necessary  condition  for  optimality 
gives  the  following  conditions. 


If  xk  >  0  then  3k.  =  0  and  X.,  =  J  Df  (y  )  and 
ij  ij  ij  fok  m 

mtPk. 


if  xkj  ■  0  then  6kj 


l  Df  (y  )  -  X  >  0  . 
**  mm  ij 


mcp 


ij 


The  above  procedure  satisfies  these  conditions.  X^  is  the  length 


of  the  shortest  chain  between  ij  . 
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The  following  theorem  establishes  the  finite  convergence  of  the 
procedure.  In  practice  a  very  quick  convergence  can  be  obtained. 


Convergence  Theorem  3.1: 


For  concave  nondecreasing  cost  functions  the  algorithm  described 
above  converges  to  an  e-optiinal  point  in  a  finite  number  of  steps. 


Proof : 


We  have  already  shown  that  if  the  above  algorithm  terminates,  it 

does  so  at  an  €-optimal  point.  Also  we  have  shown  that  each  iteration 

of  the  procedure  involves  solving  a  shortest  chain  problem  where  the 

initial  distance  matrix  has  nonnegative  elements.  Hence,  the  number  of 

3 

steps  in  each  iteration  is  finite  (-  N  where  the  number  of  nodes  is  N) . 
So,  finite  convergence  of  the  algorithm  is  achieved  if  the  number  of 
iterations  is  finite.  This  is  established  by  proving  Property  (a)  below 
for  any  nondecreasing  concave  cost  function. 


Property  (a) : 

£ 

If  the  algorithm  generates  two  consecutive  distinct  points  X  and 

£+1  £4.3  £ 

X  (X  is  obtained  by  one  iteration  starting  from  X  )  then  the 

total  cost  Z(X)  strictly  decreases,  i.e.,  Z(X^)  >  Z(X  +* )  .  This  is 

proved  as  follows. 

£  £+3 

Let  the  flow  values  in  each  arc  m  be  y^  and  y^  corresponding 

£  £+1 

to  two  consecutive  points,  X  and  X  (the  feasible  region  in  the 

£4-1  £ 

chain  space).  X  is  obtained  from  X  by  the  change  of  a  certain  set 

of  flow  chains  between  source-sink  pairs.  Let  one  such  chain  be  changed 
a  * 

from  p..  to  p,.  .  Since  p .  .  is  the  shortest  chain  using  y  flows 

rij  ij  rij  °  /ra 


in  each  arc  m  : 


56 


l  3\(>V)  i  l  B\(yl)  (I£  thr  derivative  Dfm(y‘)  exists, 
"£Pij  Btp'j 

then  the  inequality  is  also  true  when  DTf^y^j  *s  replac®d  by  )j . 

•k 

If  c  amount  of  flow  is  redirected  from  p..  to  p..  then  the 

ij  ij 

resulting  saving  is : 

«■(  i  biM)-  K  B+f»(y»)| b 0 

(“»«  nt',lj 


^because  from  concavity  D  f^y^j  >  D+f^y£j  where  =  p^ 


p .  .  0  p .  . 

ij 


*  *  \ 

and  p^  =  p^  -  p„  f)  1  .  Also  from  concavity 


_* 

’ij 


l  D+f  (y£  -  e\  >  l  D+f  ^y£)  >  £  D+f  (yl)  >  J  D+f  (y£  +  e\ 

u  m\  m  /  =  ~  m\  mf  =  L.  m\  m/  =  L1t  tn\  m  / 

'ij  mePiJ 


mCpij 


mep 


mcp .  . 

ij 


So,  a  further  saving  is  obtained  if  a  further  e  amount  is  redirected, 

and  the  maximum  saving  is  obtained  if  entire  flow  r^  is  redirected 

* 

from  p  to  p  .  A  similar  redirecting  of  other  flows  through  the 

£+1 

shortest  chains  reduces  the  total  cost  value  and  flow  values  y 

m 

£  £+1 

are  obtained,  hence  Z(X  )  >  Z(X  )  .  But  equality  is  only  achieved  when 

there  is  more  than  one  shortest  chain  and  cost  functions  are  flat  in  the 

region  involved  in  the  transfer  of  flows  or  when  flow  values  remain 

1  i 

unchanged.  In  both  cases  the  algorithm  stops  at  X  .  So,  if  X  and 

t+ 1  p  p  i 

X  are  distinct  points  then  Z(Xfc)  >  Z(Xfc  )  . 

I 

But  in  the  algorithm  described  above  the  sequence  generated  by  {X  } 
is  the  extreme  points  of  the  polyhedral  set  defined  by  (1.1b)  and  (1.1c). 
Since  the  polyhedral  set  has  a  finite  number  of  hyperplanes  it  has  a 
finite  number  of  extreme  points.  Nondecreasing  Property  (a)  guarantees 


that  no  two  elements  of  X  correspond  to  the  same  extreme  point.  Hence, 


J 
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the  number  of  elements  in  any  sequence  generated  by  the  iterations  of 
the  algorithm  is  finite.  Thus  finite  convergence  is  always  achieved. | | 

3.4  Modified  Procedures  for  Differentiable  Cost  Functions 

For  functions  having  derivatives  at  all  points  (i.e.,  differentiable), 
all  {-optimal  points  are  local  optimal  points.  In  this  case  the  following 
modified  procedures  are  sometimes  used  to  avoid  some  local  optimal  points 
and  get  near  the  global  optimum.  (i)  In  1st  modified  procedure  if  a 

£4-1  _£ 

point  X  is  obtained  starting  from  a  point  X  then  the  next  iteration 

—  04-1  _  a  oj.1 

is  started  from  a  point  X  *  0X  +  (1  -  8)X  (0  <  8  <  1)  .  However, 

Property  (a)  may  not  be  preserved.  Similar  arguments  as  in  the  proof  of 

-£  £+1 

Property  (a)  can  show  Z(X  )  >  Z(X  )  ,  but  not  necessarily 
-£  -£+l 

Z(X  )  >  Z(X  )  .  However,  if  8  is  chosen  properly  in  each  step  to 
maintain  monotonicity  then  only  convergence  Theorem  3.2  below  is  applicable, 
(ii)  In  2nd  modified  procedure  a  convex  combination  of  cost  is  used  for 
new  starting  point  ^i.e.,  here  cost  on  an  arc  m  =  B*Dfm(Xfc)  + 

(1  -  8)*Dfm(X  ) J  similar  arguments  as  in  the  proof  of  Property  (a) 

—  0  (4-1  —{4-1  _  0  94.1 

can  show  Z(X  )  >  Z(X  )  .  Now,  Z(X  )  =  8*Z(X  )  4-  (1  -  6)*Z(X  )  < 

.  B'ZfX*')  +  (1  -  BJ’ZCX*")  =  ZfX8-)  .  So,  Property  (a)  is  true  here. 

Convergence  Theorem  3.2: 

For  differentiable  concave  nondecreasing  cost  functions  the 
modified  (i.e.,  2nd  modified  procedure  and  a  special  case  of  1st  modified 
procedure  where  8  is  such  that  Property  (a)  is  satisfied)  algorithm 
converges  to  an  {-optimal  (also  local  optimal)  point. 

Proof : 


Here  termination  procedure  is  same  as  main  Algorithm  3.2.  So,  if 
the  algorithm  terminates  it  does  so  at  {-optimal  point.  Since  the  cost 


>8 


functions  are  differentiable,  € -optimal  point  is  also  a  local  optimal 

point  (by  Lemma  3.1).  The  shortest  path  procedures  in  each  iteration 

£ 

are  finite.  The  sequence  (X  }  generated  does  not  necessarily  have 

finite  number  of  elements.  However,  convergence  (not  necessarily  finite) 

2 

is  established  by  proving  the  following  Properties  (b)  and  (c)  in 
addition  to  Property  (a).  Consider  the  algorithm  to  be  a  mapping  M 
which  maps  from  a  feasible  region  of  flow  X  to  X  itself. 

Property  (b) : 

The  feasible  set  X  or,  at  least  the  subset  in  which  the  sequence 

,  i. 

IX  }  generated  by  the  algorithm  lies,  is  compact. 


Property  (c) : 


The  map  M  is  closed.  The  properties  are  proved  as  follows. 


(b)  The  feasible  region  in  X-space  is  defined  by  the  linear 

inequalities  in  (1.1b)  and  (1.1c),  hence  they  form  a  convex  set  which  is 

closed  but  which  could  be  unbounded.  In  the  solution  procedure  only  the 

required  amount  of  flows  r^  is  sent  through  the  shortest  chain.  Thus, 

£, 

the  flow  values  on  any  chain  in  a  point  X  in  the  sequence  generated  by 

the  algorithm,  is  bounded  by  r  .  =  Max  {r. . }  .  The  set  is  bounded  (since 

St  tj  i J 

*ach  element  of  X  is  bounded).  The  set  is  closed  and  bounded,  and  thus 
:ompact. 

(c)  The  closedness  property  of  map  M  will  be  proved  only  for 
unctions  with  derivatives,  which  is  the  case  here.  X  and  0  are  two 


Properties  (a),  (b)  and  (c)  are  needed  by  Zangwill's  theorem  on  algorithmic 
convergence  (Page  91,  Ref.  Z-2).  This  theorem  is  used  to  prove  the 
convergence  of  the.  map  M  defined  below. 

Closcdness  is  a  property  of  a  point  to  set  map.  Definition:  The  closed 
map  M  is  such  that  >  X°  ,  0^  e  M(X*')  ,  and  0*  -*  0°  Implies 

9  °  c  M(X°)  . 
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notations  used  to  define  a  feasible  vector  in  chain-space.  The  polyhedron 
S  defined  in  Chapter  1  for  X  is  also  defined  for  6  ,  i.e.. 


1 

Ik  I 

P  V 

y  k  k 

s  - 

0  ■  {  ij} 

)  d .  .  >  r . .  V  source-sink  pair  ij  ,0.,  >0 
k“l  ij  =  ij  K  ij  - 

If  there  are  n  source-sink  pairs  with  positive  requirements,  then  only 

n  elements  of  X  or  6  corresponding  to  an  extreme  point  are  positive. 

The  process  of  finding  the  shortest  chain  by  using  the  derivative  of  the 

£ 

cost  function  at  the  flow  value  on  each  arc  at  X  ,  corresponds  to  a 
local  search  operation  to  find  a  vector  0  which  minimizes  the  linear 

l 

approximation  of  Z(X)  in  the  vicinity  of  X  ,  i.e., 

Minimize  Z(X£)  +  DZ(X£)(0  -  X£) 

£  £ 
Subject  to  6  e  S  ,  where  DZ(X  )  is  the  derivative  of  Z(X)  at  X 

Since  X  is  constant,  the  above  problem  is  the  same  as  the  linear  program: 


Minimize  DZ(X£)-0 
Subject  to  0  e  S  . 


I  l 

If  an  extreme  point  0  solves  this  problem,  then  6  defines  a  new  set 

III 

of  shortest  chains.  The  direction  vector  d  ■  0  -  X  defines  the 

direction  of  improvement.  The  procedure  of  finding  a  point  X  of  maximum 

l  l 

improvement  of  the  cost  function  in  this  direction  is  to  find  X  *  X  +  i*d 

l  l 

which  minimizes  Z(X  +  r*d  )  subject  to  0  <  t  <  1  .  Since, 

X  *  X£  +  r*d£  ■  (1  -  t)X  +  t*0£  ,  the  vector  X  is  obtained  from  X£  by 


transfer  of  (fr^)  flow  from  the  path  corresponding  to  X1  to  the  path 

£ 

corresponding  to  0  for  various  ij  .  Because  of  the  concavity  (as 


discussed  before)  the  maximum  Improvement  occurs  when  t  ■  1  .  So, 
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14-1  111  4 

X  *  X  4-  d  =  6  and  the  map  M(X)  is  defined  as: 

M(X)  -  / 0 ;  Minimize  DZ(X)*0'  *  DZ(X)*0l  . 

I  0’eS  f 

£  q 

To  show  M(X)  is  a  closed  map  one  has  to  check:  Conditions  X  -*  X  , 

and  0*  -*•  6°  where  0*  c  M(X*)  ->  6°  e  M(X°)  .  This  is  equivalent  to 

showing  Minimum  DZ(Xo)*0'  *  DZ(X°)*0O  .  Consider  the  inequality: 

0'eS 

IDZCX1)-©'  -  DZ(X°)*8'|  <  |DZ(X*)  -  DZ<X°)l*je'|  <  6^‘c 

where  [ 0 1 |  <  c  because  of  the  compactness  of  the  feasible  set  and 
6^  -*•  0  as  X*  -»  X°  . 

Since  the  difference  is  uniformly  bounded  in  01  by  6^*c 

|Min  DZ(X*’)-0'  -  Min  DZ(X°)«0'|  <  6  -c  or, 

0’  0’  =  4 

|DZ(X*')-0*’  -  Min  DZ  (X°)  •  0 '  |  <  6  -c  ,  (because  6*  e  M(X *))  . 

0'  1 

Taking  the  limit  over  i  ,  c  being  a  finite  number 

|DZ(X°)«0°  -  Min  DZ(X°) *0 ' |  <  0  => 

O' 

DZ(X°)*0°  =  Min  DZ (X°) • 9 '  =>  0°  e  M(X°) 

e ' 

■>  M  is  a  closed  map. 

Properties  (a),  (b)  and  (c)  fulfill  all  the  conditions  of  Zangwill's 
heorem  on  algorithmic  convergence.  Here  the  solution  set  is  the  set  of 

The  map  M  is  described  over  the  continuous  variable  X  .  However,  it 
uses  only  extreme  points.  Hence,  it  goes  through  only  discrete  points 
of  the  set  and  corresponds  to  an  iteration  of  the  algorithm. 
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e-optimal  points,  so  either  the  algorithm  stops  at  a  solution  or  the 
limit  of  any  convergent  subsequence  is  a  solution. | | 

3.5  Strategies  to  get  Near  the  Global  Optimal  Point 

The  major  drawback  of  this  procedure  is  that  it  terminates  at  local 
optimal  points  and  we  have  no  idea  how  much  less  the  global  optimal 
cost  might  be.  We  have  tested  computationally  various  heuristic  approaches  to 
get  near  the  global  optimal  point.  The  different  approaches  are  as  follows. 
Approach  (a)  is  suggested  by  B  .  Yaged,  and  he  tested  it  emperically. 
Approaches  (b)  and  (c)  are  new. 

(a)  Using  a  Conv ex  Combinatio n  of  Flow  Vectors  or  Cost  Vectors 

£  £4-1 

If  X  and  X  are  extreme  flow  vectors  in  the  £th  and  (l  +  l)st 

iteration,  then  the  (SI  +  2)nd  iteration  is  started  from 
-94-9  9  9+1 

X  =  ex  +  (1  -  e)X  where  0  <  B  <  1  .  When  P  is  near  1,  the 
change  of  flow  is  small  between  two  iterations.  Sometimes  this  approach 
helps  in  climbing  out  of  the  valleys  of  the  local  optimals  to  get  at  the 
global  optimum.  There  is  no  one  B  which  gives  good  results.  The  test 
on  different  problems  shows  that  usually  the  number  of  iterations  increases 
when  such  a  convex  combination  is  used. 

Sometimes  a  combination  of  marginal  costs  is  used  instead  of  a 

£  £4-1 

combination  of  flows;  i.e.,  cost  on  arc  m  =  0*Df  (X  )  +  (1  -  0)  *Df  (X  ) 

m  m 

is  used.  Bernard  Yaged  used  different  0  values  to  get  better  local 
optituals  and  lie  reported  an  acceptable  range  of  0  in  the  interval 
(.8,  1.2).  The  use  of  this  over-relaxation  procedure  (i.e.,  0  >  1)  may 

be  advantageous  for  some  problem,  but  in  some  problems  if  a  0  >  1  is 
used,  then  some  elements  of  the  distance  matrix  may  be  negative.  So  some 
special  precaution  is  needed  to  solve  the  shortest  chain  problem. 


Sometimes  instead  of  using  marginal  cost  (i.e.,  Df  (y  )),  average 

mm 


cost  is 


is  used  ^i.  e. , 


UO' 


.  Here  the  proof  of  convergence  of  the 


algorithm  may  not  be  established.  However,  use  of  the  average  cost  value 
initially  is  sometimes  helpful  in  getting  a  good  starting  value  for  using 
marginal  cost,  particularly  when  the  cost  function  has  a  jump  at  0. 
Sometimes  a  combination  of  average  and  marginal  cost  gives  good  results. 
Each  of  the  methods  have  been  tested  computationally  for  different 


problems,  but  it  is  difficult  to  specify  which  works  where. 


(b)  Specializations  of  Step  0  of  the  Algorithm,  a  Systematic  Search 
for  Better  Starting  Points 

Different  strategies  are  employed  to  get  better  starting  points. 

(i)  Initially  a  local  optimum  is  obtained  using  the  algorithm  above  - 
then  to  get  a  new  starting  point  for  the  algorithm,  all  the  arcs  in  v.’hich 
the  flow  value  is  lower  than  a  certain  amount  are  made  very  expensive  - 
this  results  in  using  the  arcs  whose  flow  value  is  large  more  effectively 
and  not  using  arcs  with  low  flow  values.  (ii)  To  search  in  a  more  spread 
out  area,  arcs  having  positive  flow  values  for  one  local  optimal  point  are 
sometimes  blocked  to  get  a  completely  different  local  optimum. 

The  procedures  (i)  and  (ii)  are  utilized  around  the  best  local  optima] 
value  obtained  so  far,  to  determine  whether  it  is  better  to  search  in  the 
vicinity  of  the  best  local  optimal  value  (i.e.,  procedure  (i))  or  farthest 
away  from  the  best  optimal  value  (by  procedure  (ii)). 

From  the  test  results  it  is  recommended  that  procedure  (ii)  should 
be  used  initially  to  get  a  few  very  different  local  optimal  points.  Then 
procedure  (i)  should  be  used  to  search  near  the  best  local  optimal  value 


so  far  obtained. 
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(c)  Random  Starting  Point  Selection  Strategy  and  a  Sequential 

Sampling  Plan 

In  this  method  (a  similar  method  is  suggested  for  travelling  salesman 
problem  in  Ref.  [R-2])  the  global  optimal  cost  is  estimated  from  a 
series  of  observed  local  optimal  costs.  The  local  optimal  costs  are 
generated  by  selecting  random  starting  points  (not  necessarily  within 
the  feasible  region)  for  the  above  algorithm.  Suppose  that  a  sample  S 
of  size  n  of  local  optimal  points  having  values  c2>c2»  •••>  cn  has 
been  determined.  We  shall  give  a  procedure  for  estimating  the  global 
optimal  cost  given  this  information. 

We  have  no  a  priori  knowledge  of  the  distribution  of  local  optimal 
costs,  so  it  is  assumed  to  be  uniform  between  a  and  a  +  0  .  Given 
c^,  . ..,  c^  we  wish  to  estimate  a  and  0  .  (This  uniformity  assumption 
in  the  case  of  a  total  lack  of  knowledge  has  the  approval  of  both  of  the 
Bayesian  schools  [R-l] ,  i.e.,  necessarists  (Jeffreys  [J-2])  and 
subjectivists  (Savage  [E-2,  D-1J).  Jeffreys  argued  for  the  legitimacy 
of  using  a  uniform  distribution  in  case  of  a  total  lack  of  information 
by  quoting  that  "Bayes,  in  his  great  memoir,  repeatedly  says  that  the. 
principal  (i.e. ,  assigning  equal  probability  or  assumption  of  uniform 
distribution)  is  to  be  used  only  in  case  where  we  have  no  grounds  for 
choosing  between  the  alternatives."  Edwards,  Lindman  and  Savage  have 
shown  in  a  theorem  that  under  the  assumption  of  a  uniform  distribution 
of  parameters,  the  calculated  approximate  a  posteriori  distribution  agrees 
closely  with  the  actual  a  posteriori  distribution.  (Here  the  uniform 
distribution  is  assumed  of  the  parameter  rather  than  of  the  distribution 
function  itself.  However,  we  can  think  of  the  distribution  function 
itself  as  a  parameter.)  So,  the  basis  of  the  assumption  is  philosophically 
coherent  with  the  Bayesian  approach.) 
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The  estimate  of  a  is  the  estimate  of  the  global  optimal  solution. 
The  elements  of  sample  S  are  independent,  and  each  is  from  a  uniform 
density  function  between  a  and  a  +  6  .  So,  the  likelihood  function  L 
is  defined  as  the  following  function  of  the  sample  S  ,  a  and  6  : 


L(S;«,6) 


a  < 


<  a  +  9 


=  0  ,  otherwise. 


0 


So,  the  maximum  likelihood  estimates'1  of  a  and  0  are  respectively 
minimum  (c^,  ...,  cr)  ,  and  {maximum  (c^,  ...,  c^)  -  minimum  (c^ ,  ...,  cn) }  . 
Let  T  and  U  be  respectively,  minimum  and  maximum  elements  of  the 
sample  S  .  Then  the  distribution  functions  of  T  and  U  ,  F(t)  and 
G(u)  are  as  follows: 


1  -  F(t)  =  Prob  (T  >  t)  =  -  q - for  a  <  t  <  a  +  0 

density  function  f(t)  =  =  —  (a  +  9  -  t)n  ^  for  a  <  t  <  a  +  0 

dt  Qn  =  = 


=  0  ,  otherwise. 


Expected  values 


Definition:  The  maximum  likelihood  estimate  of  8  based  on  a  random  sample 

S  =  (x, ,  ....  x  )  is  that  value  of  0  which  minimizes  L(S;0) 

1  n 

(=f(x.,9)  •  ...  •  f(x  ,9)  where  f(x.,9)  is  the  density  function  of  x, 

1  n  l  l 

with  parameter  9)  considered  as  a  function  of  9  for  a  given  sample  S  . 
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E(T  -  a)  *  /  (t  -  «)f(t)  dt 


(t  -  a)  /0 

0 - nl“ 


i 

3  J ” x(l  -  x)n-1 


n6*B(2,n)  = 


n  +  1 


where  x  =  —  —  a  ,  B(2,n)  is  Beta  function  with  parameters  2  and  n  and 

_  U 

B(a,b)  =  ^ ■  ,  and  [a  =  a*  ja  -  1  .  (j  denotes  a  Gamma  function.) 

la  +  b 


G(u)  =  Prob  (U  <  u) 


for  a<u<a+0 


density  function  :  g(u)  =  “  “•  (u  -  n)n  *  for  a  <  u  <  a  +  0 

© 


=  0  ,  otherwise 

a+0  a+0 

f  (u  "  a)g(u)  du  =  f  (u  -  a)  — 

J  en 


—  / 


(u  -  a)°  ^du  = 


a+0 


E(U  -  T) 


3  dx  =  tttt  ’  (wbcrc  x  =  u~v^) 

■'o 

/  n  J  \  „  n  -  1  „  ,,/U  -  T\  e 

=  (rrr  -  ^TTT )'°  =  TTTT’0  ;  or’  H^t)  "  TTT7 


©  _  a  ,  , 

,/\Ll 

_I\ 

n  +  1 

\n  - 

l) 

,  so , 

HT  n  -  J  /  °  • 

Hence,  ^  ^ j  is  an  unbiased  estimator  of  a  .  The  variance  of 

this  estimator  gives  an  idea  of  how  good  this  estimator  is. 
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Variance  (t  -  =  e(t  -  “-j  "  E(T  “  =  E(~^l‘T  “  n  -  1 

2 

=  ( — £__}  *E(T2) - — — «*E(U*T)  +  - - - »*E(U2)  -  a2 

Vn  “  11  (n  -  l)2  (n  -  l)2 


E(T2)  =  E(T  -  a)2  +  2*E(T)«a  -  a2  =  /  (t  -  a)2  (6  +  a  -  t)n  Edt 

l  3 


„/  e  \  2 

+  2r  +  ttt r  - a 


=  n02 J' X2(l  -  x)n-1-dx  +  a2  +  *  n62*B(3,n)  +  a2  +  -| 


2  ,  29-ct 

n  +  1  * 


where 


t  -  a\ 
x  -  —r~) 


nO2!!  fn  .  2  .  29-ct _ 29? _  .  2  29-ct 

I — — —  a  n  +  1  (n  +  2)  (n  +  1)  n  +  1 

n  +  3 


a+9 

E(U2)  =  E(U  -  a)2  +  2-E(U)-a  -  a2  =  /  (u  -  ot)2-n^U  -du  + 

J  0n 

a 


+  2(“ + ^r)a  - 


2  /  n+1  2  2n6 • a  n9  .  2  2n6-u 

=  n9  lx  *dx  +  a  +  — ~  +  a  +  —  -  —  , 
/  n-tln+2  n+i 


[where 


u  -  ct\ 

X  =  — )  * 


The  joint  distribution  function  of  T  and  U  is  given  by  H(t,u) 


H(t,u)  =  Prob  (U  <  u,T  <  t)  =  Prob  (U  <  u)  -  Prob  (U  <  u,T  >  t) 

=  G(u)  if  a<u<t<ct  +  0  [since,  P(A)  =  P(A  A  BC)  + 


+  P (A  A  B)] 


=  G(u)  -  — -j  if  a<t<u<a+0. 

2 

density  function  h(t,u)  =  =0  if  t  >  u 


JL  f_  JL  (u  _  t)n_1'l  -  n(n  ~-U  f..  _  ^n"2 


r-~ 

L  en 


(u-t)  ,a<t<u<a  +  0 


=  0  ,  otherwise 


E(UT)  = 


a+0  u 

//’ 


a+0  u 


•t*h(t,u) *dt*du  = 


//*= 


-  1)  ,  ,n-2  .  ,  , 

-  (u  -  t)  u*t*dt*du 


a+0  u 


a+9  1 


j ».y  t-cu  -  t)n'2-dt-d» .  ss&jii  /  j 


a  a 


a  0 


[(x(u  -  a)  +  a)  (1  -  x)n  2(u  -  a)n  ^dx]du 


UTO 

_  n(n  -  1)  r 

on  J 


[  (u  -  a)n*u*8(2,n  -  1)  +  a(u  -  a)  *u»B(l,n  -  l)]du 


^where  x  =  ~  ,  t  =  x(u  -  a)  +  a  ,  and  u  -  t  =  (1  -  x)  (u  -  a 


=  n(n  -  1) 


•  B(2,n  -  1)  f  (-U-  —  °  j  ’u* 
a 

/u  -  a\n  ^  du 

l  e  )  ’u’  e 


du  +  a*  8(1  ,n  -  1) 


=  n(n  -  1) 


.  H KU  .  f  xn(0x  +  a)  •  Odx  +  .  f 

In  +  1  Jn  fn  -L 


n-l/G 

x  (0 


x  +  a) • dx  ^ 


where  x 


u  -  a  \ 
0  / 
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=  n(n  -  1) 


i  /  e2  a© 

\  +  «  . 

/  0  .  «\ 

1 

_n(n  -  1)  \n  +  2  n  +  1 

/  n  -  1 

\n  +  1  n/_ 

n  + 

•02  + 


a*0  +  a 


Hence, 


U  -  T) 

l  .  /  n  \ 

1  1  29 

n  -  1 ) 

‘  \n  -  lj 

1  \  (n  +  2)(n  +  1) 

,  2  26 
+  a  + 


2n 


(n  -  1) 


t+~2  °2  +  a‘6  +  “2) 


(n  -  1) 


/  n02  .  2  .  2n6-a\  2 

? v«ttt  +  a  +  vn:)  - a  = 


n0 


(n  -  1)  (n  +  1)  (n  +  2) 


Since 


E(U  -  T)2  =  E(U2  -  2UT  +  T2)  =  E(U2)  -  2E(UT)  +  E(T2) 

/  nO2  2  2n6  )  J  2  62  \  (  2@2 

=  \vn  +  a  +  rTTT'01/  ~  2\0,a  +  a  +  im) +  [jirriHim) 

2  ,  20ct  \ 

a  +VT~i) 

- n(n_ 1) _ _,02  _ 


(n  +  2)  (n  +  1) 

(U  -  T  \ 

T  -  -  _  -j-  1  given  by: 


n(n  +  2)  (n  4-  1) 


n(n  -  1) (n  -  1) (n  +  1) (n  +  2) 


(U  -  T)  ‘ 


-(Ki) 


- -  is  the  estimate  of  the  standard  error  of  the  estimator.  If  this 

n  -  1 

term  is  small,  the  variance  of  the  estimator  is  also  very  small. 

Instead  of  assuming  a  uniform  distribution  (i.e.,  6(1,1))  for  the 


local  optimals  we  could  assume  any  6  distribution.  If  we  assume 
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have  6(1,2)  distribution,  similar  calculations  can  be  made 


to  show  that  the  unbiased  estimator  of  a  (the  global  optimal  value)  is 

<T - t  .  .  _l_i  „here  c - Sins — . 


where  C 


If  the  distribution  is  assumed  to  be  6(2,1)  then  the  estimator  of 


{U  -  T  ) 

T  -  C  — x - (  . 


A  sequential  sampling  procedure  can  be  used  to  estimate  the  global 
optimal  value  where  a  trade-off  is  obtained  between  the  cost  of  computation 
and  gain  due  to  expected  reduction  in  global  optimal  value.  Assume  C 
to  be  the  cost  of  finding  one  more  local  optimum.  This  is  a  function  of 
the  number  of  iterations  (on  an  average)  to  get  to  a  local  optimum  and 
the  computing  time  per  iteration. 

The  estimates  of  a  and  a  +  6  when  a  sample  of  size  nQ  is  obtained 


are  \  T 


U  -  T  )  (  U  -  T  ) 

n  n  /  I  n  n  ( 

o  o'.  ,  /  ,,  .  o  o  ' 


n  -  1 
o 


and  )  U  + 


n  -  1 
o 


respectively  (where  U  and 

n 


T^  are  values  of  U  and  T  when  the  sample  size  is  n) .  Assume  that 
the  (n^  +  l)st  sample  will  have  a  uniform  distribution  between  the 
estimated  values  of  a  and  a  +  9  .  Then  the  expected  value  of  T  . 


is  computed  as  follows: 
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Let  p/Tn  +j\  he  the  estimated  density  function  of  . 

\  no  /  no 

1  n0  -  1 

fv1)  ‘ ?  ‘  <”<,  +  »(\  -  \)  ' 


Then 


n  +1 
o 


x  if  o  <  x  <  T 
»  =  n 


T  if  T  <x<a  +  6 
n  n  = 

o  o 


s(  V)  ■  I 


U  -T 
n  n 

U  +  — -2. 
n  n  -1 
o  o 


x*dx  + 


U  -T 
n  n 

r - P-r-°- 

n  n  -1 
o  o 


/ 


T  *dx 
n 

o 


■  T 


U  -  T 
n  n 
o  o 


n  2(n  +  l)(n  -  1)  * 
o  o  o 


Therefore,  given  a  sample  of  size  nQ  ,  one  of  the  following  plans  can 
be  used. 


(i)  Determine  a  function  G(T)  ,  monotonically  increasing  in  T  , 
which  gives  a  measure  of  cost  when  the  so  far  obtained  minimum  local 
optimum  value  is  T  .  We  would  like  to  have  the  sum  of  the  cost  of 
computation  and  G(T)  be  as  small  as  possible.  Now  the  estimated 
reduction  in  G(T)  when  the  (nQ  +  l)st  local  optimum  is  found  is 
|g(T)  -  G  (E(V))]  and  the  cost  of  taking  one  more  sample  is  C  . 
Therefore  the  plan  should  be:  If 


!°(\)  -  c(\  - 


<  C  , 


no  further  sample  is  taken.  The  best  local  optimum  so  far  obtained  is 


/  U  -  T  \ 

I  n  n  V 

accepted.  I  Note  that:  E/T  =  T  -  -57 - '•>'<'/ — ~ — • 

\  IV1/  no  2(no  +  D(no  -  1)/ 


’(\)  ’  G(\  ’  2(no  +  1)<no  "  1}){ 


If 


>  C  , 


then  another  sample  is  taken  and  the  criterion  function  is  recomputed. 

The  maximum  number  of  samples  N  taken  is  bounded  because  (U  -  T  ) 

n  n 


is  bounded  by  0  .  Hence, 


U  -  T 
N  N 


becomes  small  for  large  N 


2(N  +  1) (N  -  1) 

making  the  difference  between  values  of  G(*)  at  T  and  E 


n  (Tn  +l) 

o  \  o  / 


small  enough  to  be  less  than  C  . 

(ii)  Define  a  loss  function  L(e)  ,  a  monotone  increasing  function 
of  e  ,  which  gives  the  cost  incurred  when  the  estimated  expected 
decrease  in  the  global  optimal  value  by  taking  one  more  sample 
^i.e.,  T^  -  E^T^  +lj^  is  e  .  In  this  case  the  plan  should  be:  If 


/  U  -  T  \ 


no  further  sample  is  taken.  The  best  local  optimum  so  far  obtained  is 


accepted.  ^Notc  that  [t^  -  E^T^J  "  T^VjXn^-  1)  J  ’  lf 


U  -  T 
n  n 


/  U  -  T  v 

( _ _ ) ,  c 

\2<"o  +  ‘><"o  -  »)-  ’ 


then  another  sample  is  taken  and  the  criterion  function  is  recomputed. 
The  maximum  number  of  samples  N  is  bounded,  because  if  the  procedure 


stops  at  N  then. 
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/  un~tn  \  .  t/un-i  -  Vl\ 

J\2(N  +  1)(N  -  1)/  L  3nd  L\2(N)(N  -,2)/ 


>  C 


So, 


N(N  -  2)  < 


U  -  T 
N-l  XN-1 

2L"1(C) 


2L~1(C) 


(Since  L(*)  is  a  monotone  increasing  function  and  its  inverse  L  *"(•) 
exists.)  Or, 


(where  6 
is  bounded. 


N(N  -  2) 


< 


n  +  1 
o 

n-l 


o 


U  -  T 
n  n 
o  ' 

L_1(C) 


is  replaced  by  its  estimator)  N  is  bounded  because  N(N  -  2) 
Hence,  this  procedure  stops  after  a  finite  number  of  samplings. 


This  sequential  sampling  approach  has  been  tested  using  C  and  L 
as, respectively,  an  increasing  linear  functions  of  t  ,  the  average 
number  of  iterations,  and  e  ,  the  estimated  expected  decrease  in  the 
value  of  T  by  one  more  sample  (i.e.,  C  =  at  and  L  =  be  where  a  and 
b  are  specified  coefficients) .  In  the  solution  procedure  the  value  of 
a/b  used  is  between  30  and  40.  A  comparison  was  made  of  this  sequential 
procedure  with  the  global  optimal  search  procedure  (of  Chapter  4)  by 
determining  the  computation  time  of  the  global  search  procedure  and  actual 
error  and  computation  time  of  this  sequential  search  procedure.  About  25 
problems  of  varying  sizes  (number  of  nodes  between  6  and  35)  have  been 
solved  using  both  the  global  search  procedure  of  Chapter  4  (with 
approximately  2Z  error  value)  and  the  sequential  search  procedure  of  this 
chapter  (using  a/b  *  35).  The  estimated  error  upon  termination  of  the 


sequential  search  procedure  was  on  the  average  around  3%  while  the  actual 
error  was  about  8  -  10%  (casting  some  doubt  upon  our  use  of  a  uniform 
distribution  for  the  values  of  local  minima).  The  computation  time  for 
the  sequential  search  procedure  was  around  50  -  60%  that  of  the  global 
search  procedure  of  Chapter  4.  So,  computationally,  the  sequential 
search  procedure  seems  to  be  attractive.  However,  we  cannot  compare 
these  two  procedures  for  problems  with  a  greater  number  of  nodes  because 
of  limitations  on  the  global  search  procedure  caused  by  restricted 
computer  memory  size  and  computational  speed. 


74 


CHAPTER  4 

In  this  chapter,  a  solution  procedure  will  be  discussed  which 
converges  to  a  value  as  close  as  desired  to  the  global  optimum.  This  is 
a  branch  and  bound  method.  The  procedure  is  an  adaptation  of  a  general 
procedure  due  to  Walkup  [W-l]  and  Falk  and  Soland  [F-l]  to  our  special 
problem.  The  general  procedure  requires  the  solution  of  certain  sub¬ 
problems.  These  subproblems,  for  our  problem,  have  a  special  structure 
and  we  have  developed  special  procedures  to  exploit  this.  Each  subproblem 
involves  finding  the  shortest  chains  between  all  pairs  of  nodes.  The  first 
major  difficulty  faced  in  solving  a  large  problem  by  this  method  is  the 
limited  size  of  a  computer  memory.  Some  ways  to  get  around  this  problem 
will  be  described  in  this  chapter.  The  second  major  difficulty  arises  in 
determining  the  upper  bounds  on  the  flow  in  each  arc.  An  inefficient 
upper  bound  is  the  sum  of  all  the  flows.  Different  types  of  heuristic 
methods  which  yield  more  efficient  bounds  will  be  discussed.  However ,  the 
bounds  are  very  critical,  and  choice  of  the  wrong  bounds  may  lead  to 
convergence  to  a  nonglobal  optimal  solution. 

4.1  Reformulation  of  the  Problem 

The  problem  can  be  reformulated  with  the  superfluous  extra  constraint 

(4.1c)  shown  below.  This  constraint  provides  us  with  a  range  in  which  the 

concave  function  could  be  approximated  to  develop  a  solution  procedure. 

This  range  1°  defined  below  is  critical  for  the  effectiveness  of  this 
m 

procedure. 

M 

(4.1a)  Minimize  Z(Y)  =  Y  f  (y  ) 

** ,  n:  m 
m=i 

(4.1b)  Subject  to  Y  r.  S 

M 

Y  e  n  1° 
m*=l  m 


(4.1c) 
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where 


S 


Y 


<y> 

D3 


y  =  I 


A.  .X.  . 
3  3 


13 


k-1 


x .  . 
ij 


>  r.  . 
=  iJ 


and 


and  interval 


m 


0 


r .  . 

ij 

ij  3  a  chain  between 

ij  contains 

the  arc  m  and  i>j 


At  this  point,  the  following  definitions  are  useful. 


4.2  Definitions 

The  linear  approximation  of  a  function  f  (,*)  in  the  interval 

I  =  [1  ,u  ]  is  the  function  f^inr(.)  given  by  flinr(y  )  = 
m  m  m  m  m  ni 

a  +  b  (y  -  1  )  where  a  =  f  (1  )  and  b  =  (f  (u  )  -  f  (3  ) ) / (u  -1). 

mmmm  mmm  m  mm  mm  mm 

This  is  drawn  in  Figure  4.1.  For  a  concave  function 

flinr(y  )  <  f  (y  )  in  the  interval  1  and  f^1Ur(y  )  >  f  (y  )  outside 

m  m  =  m  m  m  m  m  =  m  m 

the  interval  I 

m 

Linearization  Error 

The  difference  e  (y  )  =  ]f  (y  )  -  f^nr(y  )|  is  called  the 
mm  [_rnm  m  '  m  J 

linearization  error  of  f  (y  )  at  y  .  This  error  is  nonnegative  in 

mmm 

the  interval  I  and  nonposilive  outside  the  interval  I 
m  m 

4.3  Linearized  Version  of  the  Problem 

In  this  procedure,  instead  of  solving  (4.1),  a  linearized  version 


of  this  formulation,  given  in  (4.2),  is  solved  initially. 
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M  1  . 

(4.2a)  Minimize  Z * (Y )  =  J  f  lnr(y  ) 

,  m  m 


(4.2b)  Subject  to  Y  e  S  where  S  is  as  defined  in  (4.1) 

M 

(4.2c)  Ye  11  1° 


,  ,linr,  s  .  , 

where  f  (y  )  is  the 
m  m 

f  (y  )  in  the  interval  1° 
m  m  ra 


Let  Y 


linr 


linear  approximation  of  the  function 


.opt 


be  a  solution  of  (4.2)  and  Y  be  a  solution  of 
>linr  „.,linr.  „mixr  „  „,linr. 


//  ™  ,linr  ,  .  lmr.  mixr  ,  lmr.  ,  „opt  -/..opt. 

(4.1).  Define  Z  =  Z' (Y  )  ,  Z  =  Z(Y  )  and  Z  v  -  Z (Y  1  )  . 

Then  the  following  inequality  is  obtained 


(4.1) 


zlinr  <  zopt  <  zmixr  =  zlinr  +  E 


where  E  =  V  e  (y^inr)  and  e  /y^^nr)'s  are  linear ixali 
L .  mV  rn  /  mV  m  / 


on  errors  a  I 


m=l 


,linr 


f  lmr)  ,  ,  ,o 

<y  >  for  intervals  1 

(  m  J  m 


The  first  inequality  is  obtained 


because  within  1°  ,  f^^nr(y  )  <  f  (y  )  and  the  answer  is  within  J°  . 

mm  m  =  m  m  m 

The  second  inequality  is  obtained  because  y^nr  fs  any  feasible  solution 
of  (4.1)  and  7°^^  is  the  optimal  solution.  The  final  equality  is 
obtained  because  the  extra  amount  K  is  the  suit,  of  all  linearization 


errors . 


If  E  is  small,  then  Z  and  Z  Xr  will  be  a  good  approximation 

of  the  optimum  value.  The  maximum  value  of  E  depends  on  the  intervals 

l°'n  ;  partitioning  of  those  intervals  into  shorter  ones  reduces  the 

maximum  possible  value  of  f  .  The  general  procedure  involves  solving. 

subprobl ems  of  the  form  (4.7a),  (4.21>)  subject  to  the  additional  eonstiainl 

that  the  flow  in  each  arc  m  is  restricted  to  lie  in  a  .specified  interval 

[1  ,ti  ]  .  At  each  iteration,  two  raw  subproblems  are  generated  from  a 
mm 
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given  subproblem  by  dividing  the  admissible  interval  for  one  arc  m' 

into  a  left  part  [1  .  ,y  .]  and  right  part  [y  , ,u  ,]  and  leaving  the 

mm  mm 

admissible  intervals  for  all  other  arc  flows  unchanged.  Hence,  each 
iteration  adds  1  to  the' size  of  the  set  T  of  subproblems  where  T 

has  the  property  that  any  (admissible  by  (4.1fc))  flow  pattern  is 

admissible  for  at  least  one  subproblem.  (If  the  flow  happens  to  correspond 
to  a  point  of  division  of  any  one  arc,  it  will  be  admissible  for  two  sub¬ 
problems.)  At  iteration  p  ,  the  original  problem  has  been  partitioned 
into  p  subproblems.  To  be  able  to  uniquely  identify  a  subproblem  by 
number,  we  number  the  subproblems  as  follows.  Assume  the  subproblems  at 
iteration  p  are  uniquely  numbered  from  1  to  p  (such  is  the  case  at 
iteration  1)  .  If  now  the  subproblera  k  is  partitioned  into  two  sub¬ 
problems  by  dividing  arc  m1  at  y  ,  the  subproblem  associated  with  the 
left  hand  part  Um,.jr]  Is  numbered  k  and  the  subproblem  associated 
with  the  right  hand  part  is  numbered  p  +  1  .  Let  1^*^  be 

the  admissible  interval  of  arc  m  for  subproblem  k  at  iteration  p  . 

M  p  M 

Then  n  1°  =  U  n  I  for  all  p  .  If  we  say  that  the  kth 

i  m  ii  i  m 

m=l  k=l  m=l 

subproblem  is  solved  at  iteration  p  ,  we  mean  the  problem 


('■4a) 

(4.4b) 

(4.4c) 


Minimize 


Subject  to 


M 


Z’(Y)  =  l  fU"r  (y  ) 
m,k,p  -'m 

m- 1 


Y  c  S  ,  S  defined  above 

M 


y  c  n  i 

m=l 


k,p 


m 


where  f^J11  (y  )  is  the  linear  approximation  to 

m,k,p  m 


corresponding  to  I 


k,p 


m 


Bounding  refers  to  obtaining  the  lower  bound 

(4.4);  i.e.,  if  Yj^nr  is  the  solution  for  (4.4) 
k. ,  p 


the  cost  on  arc  m 

for  each  subproblem 
then 
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linr 


Z,(Y?':^nr)  is  the  lower  bound.  Define  Z™*Xr  =  zfy^nr\  and 

k,p  \  k,p  /  k,p  \  k,p  / 

^nr  =  Minimum  <Z,  lnri  when  p  is  the  number  of  subproblems  solved. 
s,P  k=1  2  _  l  k,p  / 


Then  the  following  inequalities  are  obtained 


(4.3) 


,Hnr  <  zoPt  <  M.niinum  (zmixr)  <  zmi 
s’p  =  =  k=l,...,p  l  k*P  >  =  s> 


mixr  „linr  ,  „ 
=  Z  +  h 


s,p 


s,p 


,  Y  /  lmr  \  ,  linr 

where  E  =  )  e  { y  }  and  y 

s,P  m\  m,s,p/  m,s  ,j 

m-l 


vector  Y 


linr 

s,p 


is  the  mth  element  of  the 


.linr 


The  first  inequality  is  obtained  because  each  Z  gives  the 

M  .  k,P 
k.  d 

lower  bound  for  the  problem  in  the  interval  II  I  and  the  minimum 


m=l 


m 


gives  the  lower  bound  for  the  entire  interval,  hence  the  lower  bound  for 

Z°P^  .  The  second  inequality  is  obtained  because  each  z!"*xr  is  a 

k,p 

feasible  solution  of  the  main  problem  (4.1),  hence  Z°pt  is  less  than  or 

equal  to  it.  The  third  inequality  is  obtained  because  2mixr  is  one  Df 

s  ,p 

the  elements  over  which  minimization  is  done,  and  the  final  equality  is 


obtained  because  E  has  taken  care  of  the  linearization  error.  In 

s,p 

the  branch  and  bound  process,  the  maximum  value  of  E  is  made  small 

s,p 

,  ,  ,  ..opt  .  „linr  .  , 

so  that  a  good  estimate  of  Z  is  Z  .  This  above  error  can  be 

s,p 

made  arbitrarily  small  in  a  finite  number  of  steps  as  shown  below. 
Moreover,  the  following  propositions  make  the  branch  and  bound  procedure 
computationally  attractive. 


Propos i t ion  4.1 : 

:l  may  not  be  necessary  to  break  all  the  rectangles  of  a  decomposition 

down  until  t'c  maximum  of  the  error  estimate  K,  for  each  one  is  as 

k,p 

smal '  as  the  desired  aceurury  of  an  answer,  say  <  . 
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Proof : 


If  any  feasible  solution  Y°  with  the  value  Z°  =  Z(Y°)  is  known 
for  the  original  problem  (4.1),  and  if  the  solution  zf'*'nr  of  (4.4)  is 

k,p 

such  that  Z  lnr  >  Z°  ,  then  the  optimal  solution  cannot  lie  in  the 

K  ,  p 


M 


rk,p 


rectangle  II  I  *  (because  its  lower  bound  has  a  greater  value  than 
m 


m=l 


the  value  obtained  at  some  other  point).  So,  r.o  further  refinement  of 


that  subproblem  is  necessary.  In  the  procedure,  the  minimum  Z™lxr  so 

k ,  p 

far  obtained  is  taken  as  Z°  and  at  each  step  all  the  intervals  for 

which  zjinr  >  Z°  -  €  are  dropped  from  further  consideration.  This 
k,P  — 

procedure  of  rejecting  some  intervals  altogether  hastens  the  process  of 
convergence. I | 


Proposition  4.2: 

The  procedure  of  the  algorithm  described  below  remains  unchanged  if 
the  restriction  (4.4c)  is  dropped;  i.e,,  the  kth  subproblem  solved  at  the 
pth  iteration  is 


(4.4a) 

(4.4b) 


M  . 

Minimize  Z'(Y)  =  V  f  ,r  (y  ) 

t  m,k,p  m 
m-i 


Subject  to  Y  e  S 


Proof : 

Let  y^n  be  the  solution  for  this  unrestricted  problem  and  define 

n  L  in  ,  ,  lin.  .  mix  _«liii.  .  . , .  .  . 

Z,  -  Z  (Y  )  and  Z  =  ZO  )  .  (Notice  the  difterencc  between 

k,p  k,p 

vlinr  ,linr  ,  „mixr  .  „]in  _lin  „mix  , 

Y  ,  Z^  and  Z^  and  Y  ,  Z^  ,  7.  .  The  last  r  Is 

dropped  for  the  unrestricted  problem.)  Here  Y^n  no  longer  lies  within 

the  interval  so  inequality  (4.5)  is  not  obvious.  However,  in  the 

M  , 

restricted  problem  the  constraint  is  Y  e  S  f)  n  1  and  in  tin- 


m=l 


unrestricted  problem  the  constraint  is  Y  e  S  . 


Since 
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S  0  II  I  ,p  C  S  ,  it  is  obvious  that. 

_i  m 

m=l 


Minimum 

M  1 
Y  e  S  fl  II  I 


Z'(Y)  >  Minimum  Z' (Y) 
YeS 


The  better  (smaller)  minimum  is  obtained  when  the  search  is  done  on  a 


larger  set  than  on  one  of  its  subsets.  So,  if  V. 


Minimum  <zj**ni 


where  p  is  the  number  of  subprobloms  solved,  then  the  above  inequality 
plus  the  first  inequality  of  the  chain  (4.5)  establishes  the  first  of 
the  following  chain  of  inequalities: 


(4.6) 


in  opt  „  ....  (mix)  _mix  .1  i 

<Z  <  Minimum  <Z,  >  <  Z  =  / 

’ » P  ~  s_i  n  l  * ,  p  J  S»P  3 » 


lln  +  E 
s,p  s,p 


Since  all  Z  are  feasible,  the  second  and  third  inequalities 

k ,  p 

follow  as  before.  The  final  equality  is  also  valid  because  E  takes 

s,p 

care  of  the  linearization  error.  Nov;  the  maximum  value  of  E  is  the 

s,p 

maximum  value  that  the  linearization  error  can  have  within 
M  ,  .  .  . 

II  1  .  Because  m  this  ease  some  y  (nit  1  i  element  of  \  )  may 


not  be  within  the  interval,  the  errors  corresponding  to  those  arc  negative. 

And  for  the  ones  which  are  within  the  interval  the  error  is  positive  and 

is  bounded  by  the  maximum  linearization  error  of  the  interval.  So,  the 

sum  E  is  bounded  bv  the  sum  of  the  maximum  linearization  error  within 

s,p 

the  intervals.  The  proposition  4.1  is  also  valid  here  because  Z. ln 

k ,  p 

does  indicate  at  least  a  lower  bound  of  the  objective  function  on  the 

M  . 

interval  II  I  .  So,  both  the  inequality  chain  (4.6)  and 

-i  m 

m=  1 

Proposition  4.1,  which  are  necessary  for  the  algorithm  to  work,  are 
satisfied  in  this  case .11 


If  the  restriction  of  bounds  on  flow  in  each  arc  is  dropped,  then 


the  solution  of  the  unrestricted  subproblem  is  accomplished  by  solving 

for  the  shortest  chain  between  each  pair  of  sources  and  sinks  with  the 

length  of  each  arc  equal  to  the  slope  of  the  linear  cost  curve  and  by 

passing  the  entire  required  flow  through  the  shortest  chain.  The  order 

3 

of  calculations  in  each  step  is  N  /2  where  N  is  the  number  of  nodes. 


to 

rik’p,rl 

and 

r — 
e- 

J 

• 

u± 

L  m  J 

L  m  J 

where  lk,p  <  r  <  uk,P  ,  then  because  of  the 
_  m  m 

concavity  of  the  cost  function  the  slope  of  the  linear  curve  for  £lk,P,rJ 
increases  and  the  slope  of  the  linear  curve  for  £r,uk,PJ  decreases  as 
shown  in  Figure  4.2.  The  two  subproblems  obtained  by  dividing  the  arc  are 
as  follows: 


(1)  Solution  of  a  problem  where  one  arc  length  has  decreased.  Here 

the  modified  solution  procedure  described  in  the  algorithm  can 

2 

be  used  to  solve  the  problem  in  -  N  steps.  Thus,  considerable 
computation  time  is  saved. 

(2)  Solution  of  a  network  problem  where  the  length  of  one  arc  has 
increased.  Usually,  this  requires  solving  the  problem  all  over 
again  in  which  case  the  calculations  are  ~  N^/2  .  However,  if 
the  same  arc  is  divided  again  and  again,  then  use  a  very  large 
length  for  the  particular  arc  and  initially  calculate  the 
shortest  chain  matrix  and  store  it.  Even  though  the  arc  length 
has  increased  in  a  particular  subdivision,  it  has  effectively 
decreased  from  the  arc  length  used  in  calculating  the  initially 
stored  matrix.  So,  the  modified  procedure  described  in  the 
algorithm,  using  the  stored  matrix  can  solve  this  new  subproblem 


Cost 
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The  algorithm  is  described  in  full  detail  below. 


4.4  Algorithm 

Step  1:  (Initialization) 

For  each  arc  m  ,  an  interval  1°  =  [0,u  3  is  determined.  (The 

m  m 

determination  of  u  ,  the  upper  bound  on  possible  flow  on  arc  m  in  the 

m 

optimal  solution  is  difficult,  and  the  convergence  of  the  procedure 

depends  on  how  well  it  is  obtained.)  The  linearized  function  Z' (Y)  for 

this  set  of  intervals  is  formed  and  the  unrestricted  problem  given  by 

(4.4a)  and  (4,4b)  is  solved  to  obtain  Y^in  ,  Z^in  (i.e..  Minimum  of 

Z'(Y)  without  the  upper  and  lower  bound  restriction  on  variables)  and 

^mix  ^  (The  subscript  indicates  that  this  is  the  1st  subproblem.)  The 

solution  procedure  requires  solving  the  shortest  path  problem  using  Floyd's 

algorithm  between  all  pairs  of  points  and  passing  the  entire  flow  through 

the  shortest  path  between  each  source-sink  pair.  Y^in  and  Z^ix  is  set 
bgg  £  best 

as  Y  and  Z  ,  the  best  solution  found  so  far.  The  information 

describing  the  intervals  1°  and  Z^in  is  stored  as  the  first  and  only 
item  of  a  list.  Set  p  ,  which  indicates  the  iteration  number,  equal  to 

1  . 


Step  2: 


Every  item  of  the  list  for  which  the  associated  value  of  Z^in  is 
greater  than  or  equal  to  (Z^eSt  -  ()  (where  €  is  the  permissible  error) 
is  removed  and  discarded.  If  the  list  is  now  empty,  go  to  Step  6. 
Otherwise,  select  the  item  from  the  list  corresponding  to  the  minimum 
value  of  Z^°  .  Assume  this  is  the  value  found  when  the  kth  subprobJem 
is  solved.  From  among  the  intervals  described  in  the  selected  item, 
choose  the  interval  I  for  which  the  error  e 


at  the  solution  point  is 
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maximum,  i.e. ,  e 


=  Max  /f  (y  ln)  -  f11"  (ylin)l  and  this  interval 
■  l  m\  m  /  m,k,p\  m  / J 


I  =  [1  ,u  ]  is  broken  into  two  intervals  1^  =  fl  ,y^^n|  and 
s  1  s’  s1  s  L  s  s  J 

1  =  y  ,u  I  .  (The  information  about  which  interval  should  be  divided 

S  S  SI 

for  each  item  can  be  obtained  in  the  same  step  when  the  Z^in  is 

evaluated,  so  that  y^in  values  for  all  arcs  need  not  be  stored.) 

m 


Step  3: 


Compute  the  linear  approximation  to  the  arc  s  based  on  the 

interval  1^  and  leave  all  other  intervals  as  in  the  kth  subproblem. 

(4.4a,b)  is  solved  to  obtain  Y^n  ,  and  Zm*X  .  Call  this 

subproblem  (k,p  +  1)  .  If  is  ^ess  than  Z^CSt  ,  Y^eSt  and 

„best  ,  .  ,  „lin  .  „mix  „ .  ,, 

Z  are  replaced  by  and  Z^  j  .  finally,  the  description 

of  the  intervals  for  subproblcm  (k,p  +  1)  and  the  value  of  Z?^n  ,  is 

k ,  p+ 1 

added  as  a  new  item  in  the  list.  (This  step  is  an  all-pairs;  shortest 

3 

path  problem  and  requires  -  N  /3  additions  and  comparisons.)  (All 
other  subproblems  ( (p  -  1)  at  most)  except  the  kth  one  at  this  stage 
will  change  its  second  subscript  to  p  +  1  from  p  because  these  sub- 
problems  will  be  used  in  the  next  iteration  only.) 


Step  4: 


2  1 

Step  3  is  repeated  using  1^  in  place  of  I  .  Call  this  subproblem 

(p  +  1  ,  p  +  1 )  .  However,  here  the  length  of  arc.  s  (suppose  it 

* 


connects  n  and  1)  is  shortened  to  d 


(This  is  b  as  defined 
m 


I 

m 

.  Here  m  =  s 

and 

I  =  l2  .) 

m  s 

The 

new 

two 

nodes  ij  is 

* 

"u 

■  V 

,  d .  +  d  + 

m  nl 

Vil 

where  d  . 

i  .1 

i  s 

the  shortest 

and 

j  when  the  arc  is 

not  slinvt  i  m  J 

•  ( 

liie  op t  i ma  1 

solution  is  obtained  with  on  1  v  N  c.ih  ubi!  ions.  ) 
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Step  5: 


Delete  and  the  description  of  the  Intervals  used  In  the  kth 

subproblem  from  the  list.  Add  1  to  p  and  go  to  Step  2. 


Step  6: 


Stop.  The  vector  Yc 


is  the  feasible  solution  to  the  original 


best 

problem  and  its  value  Z  is  within  €  of  the  optimal  value. 


4.5  Remarks 


Various  options  other  than  those  given  in  the  algorithm  can  be  used. 
The  choice  of  the  arc  to  be  considered  for  further  subdivision  can  be 
arbitrary.  Using  the  same  arc  again  and  again  for  a  few  iterations  reduces 
the  storage  requirements  and  can  be  done  efficiently.  Also,  the  sub¬ 
division  point,  given  an  interval  and  an  arc,  can  be  the  mid-point  of  the 


interval,  the  point  x  where  the  linearization  error 


{frn(x)  - 

lm  m,K,p  | 


is  maximum  or  a  point  which  divides  the  interval  such  that  the  maximum 

linearization  error  in  each  part  is  equal.  The  relative  merits  of  these 

can  only  be  tested  empirically.  However,  the  use  of  y^*n  as  the  point 

m 

of  partition  yields  good  results. 

This  branch  and  bound  algorithm  is  equivalent  to  generating  a  tree 
as  follows:  Step  1  produces  node  a  and  puts  information  about  the  node 
an  a  list.  Step  2  finds  any  pendant  nodes  of  the  current  tree  that  need 
lot  be  considered  further,  flags  them,  and  deletes  information  about  them 
from  the  list.  If  the  list  is  empty,  it  stops.  If  not,  it  chooses  a 
sendant,  unflagged  node  (call  it  node  x)  and  creates  two  descendants  of 
chat  node.  Steps  3  and  4  perform  calculations  on  the  two  new  nodes  and 
idd  information  about  these  two  nodes  to  the  list.  Step  5  deletes 
Information  about  node  x  from  the  list,  flags  node  x  and  returns  to 
'>tep  2.  The  algorithm  terminates  at  Step  6  when  all  the  nodes  are  flagged. 
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At  any  stage  of  the  algorithm,  all  the  intervals  over  which  linearization 
is  made  and  which  are  associated  with  the  terminal  nodes  of  the  tree, 
cover  the  entire  set  of  initial  intervals  1°  and  hence  the  polyhedron 
S  . 


Convergence  Theorem  4.1; 

Given  any  arbitrarily  small  value  of  €  >  0  and  assuming  that  each 

function  f  (y  )  is  concave  and  nondecreasing , ^  the  algorithm  described 
m  m 

above  terminates  in  a  finite  number  of  steps  and  correctly  produces  a 
feasible  solution  with  the  criterion  value  within  e  of  the  global 
optimal  value. 


Proof: 

Solution  of  two  subproblems  is  needed  in  each  iteration  of  the 
algorithm.  Each  of  these  subproblems  is  a  shortest  chain  problem  with  all 
nonnegative  entries  in  the  initial  distance  matrix.  Hence,  each  can  he 

3 

solved  in  a  finite  number  of  steps  (i.e.,  maximum  -  N  steps,  where 
N  =  number  of  nodes).  So,  we  will  only  have  to  show  that  number  of 
iterations  p  needed  before  we  go  to  Step  6  is  finite  for  a  given  €  , 
and  also  that  once  we  go  to  Step  6  we  have  found  a  feasible  solution  with 
criterion  value  within  €  of  the  global  optimal  value. 

Suppose  in  some  iteration  p  we  have  reached  Step  6,  i.e.,  all  the 
pendant  nodes  of  the  tree  (all  p  subproblems)  are  flagged.  According  to 
flagging  rule  Z^est  -  €  <  for  all  values  of  k  (i.e., 

k  "  1 ,  . . . ,  p)  .  Since  all  values  of  k  cover  the  entire  feasible  region 

^ Nondecrcus i ng  pr».  petty  Is  inn-ess. iry  bote  because  it  est  .i!>  I  I -.lies  the 
continuity  of  tin-  function  at  tin  right  ham!  boundary,  whi»h  is  need'll  in 
the  proof.  Also,  this  prepeity  )  uni  ant  ees  tii.it  all  tin-  entries  oi  the 
distance  matrix  used  to  solve  tin  shortest  path  problem  are  nonncgal i or . 


8 


<  M  p  M  \ 

i.e.,  II  1°  »  U  n  I  ’^1  of  problem  (4.1) 

(  m-1  k=l  m=l  / 

(“t"”  HTp))  - 


and  Z°pt  is  the  global 


iptimal  solution,  then 


Z°pt  ,  and 


best 


-  €  <  /Minimum 
*  \  k 


{*£})  i  z°Pt 


,  so  we  have  obtained  Z 


best 


within 


e  of  the  global  optimal  solution. 

Therefore  all  we  will  have  to  show  is  that  we  go  to  Step  6  in  a 

inite  number  of  iterations  p  for  any  given  €  >  0  .  Observe  that  the 

unctions  are  concave  and  nondecreasing,  hence  they  are  continuous  except 

ierhaps  at  zero  flow  value  (where  there  can  be  a  jump).  So,  except  for 

ntervals  in  the  neighborhood  of  zero  (which  we  consider  later)  the 

inearization  error  e  (•)  for  any  interval  [y  ,y  +  I  )  is  bounded  by 

in  m  m  m 

he  following  inequality: 


e  (*) 
m 


(f  (y  +  I  )  - 
mm  m 


m  m 


(d+£  (y  )-I  )  , 

v  mm  m/ 


here  D  denotes  the  right  hand  derivative.  The  first  inequality  is  due 

o  the  fact  that  f  (•)  is  nondecreasing,  and  the  second  inequality  is 

ue  to  the  fact  that  f  (•)  is  a  concave  function.  The  right  hand 

m 

»rivative  D+f  (•)  is  finite  except  at  y  =  0  (shown  in  Appendix  D) . 
m  m 

>w  at  each  iteration  one  of  the  intervals  is  partitioned  into  two,  so 

>me  I  decreases.  If  we  can  show  that  by  repeated  partitions  of  an 
m 

>terval,  I  -*■  0  ,  then  e  (•)  will  become  less  than  C/M  in  a  finite 
m  m 

imber  of  steps.  Consequently,  the  total  error  satisfies 


M 

l 

m»l 


e  (•) 
in 


<  €  at  some  finite  number  of  iteration  p  for  all  k 


nee  from  (4.6) 


zlln  <  Z°P 

s.p  - 


1  <  (  Minimum 

"  Wi . .  t  k’p>/ 


-best 
Z  < 


+  E  and  since  total  error  E 
>P  s,p  s,p  s,p 

Cst  <  Zlin  +  <  or,  ZbeSt  -  €  <  Zlin  <  Zyn 
=  s.p  -  s,p  =  k , p 


is  bounded  by  l 
V  k  =  1,  ....  p  because 
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is  the  minimum  element  of  all  the  * n  .  So,  in  finite  number  of 
s,p  k,p 

iterations  p  ,  all  the  pendant  nodes  satisfy  the  flagging  criteria.  Hence, 

we  reach  Step  6  in  a  finite  number  of  steps.  We  now  consider  intervals 

containing  0  ,  then  to  complete  the  proof  we  show  that  1  ->  0  . 

m 

The  bound  e  (*)  <  D+f  (y  )  .  1  is  not  vaiid  for  an  interval  I 

m  =  m  m  m  m 

with  left  hand  endpoint  0  if  D+f  (0)  does  not  exist.  it  appears  that 

m 

we  might  find  ourselves  reducing  the  length  of  1  without  reducing  the 

linearization  error  at  the  point  of  division  y  below  the  amount  of 

m 

discontinuity  at  0  .  This  could  happen  if  y  -*■  0  .  But  the  procedure 

m 

guarantees  the  flow  value  on  each  arc  is  either  zero  or  at  least 

r  =  Minimum  { r . . }  .  So,  such  points  do  not  exist. 

S  ij  3r . . >0 
ij 

The  only  other  possibility  Jor  the  maximum  linearization  error  to 
remain  large  for  some  interval  1^  in  subsequent  partitions  is  if  the 

partitions  are  always  very  near  one  of  the  boundary  points  (i.e.,  one 

I  decreases  very  little  in  each  iteration  and  approaches  a  nonzero 
limit).  But,  in  _his  case,  because  of  the  continuity  of  the  cost  curve, 

the  actual  error  eg(y  )  at  point  wherr  the  interval  is  partitioned 

eventually  becomes  <  €/M  .  (Because  close  to  the  boundary  point  where 
the  linearization  error  is  zero,  linear  approximation  and  actual  curve 
are  close  together.)  Since  partition  is  made  on  an  arc  s  such  that 


c  (y  )  =  Maximum  (e  (y  ) }  ,  the  total  actual  error 
s  s  m  m 

m 


1 


e  (y  )  <  M*o  (y  ) I  is  <  €  whore  e  (y  )  <  e /M  .  The  values 
mm  —  s  s  /  ss 


of  zf*  and  for  this  subproblem  is  very  close,  so  they  will  be 

a  ,  cc  l  i  ,bcst  ,  _lin  ,  ..best  „mix 

flagged  off  by  the  rule  Z  -  e  <  7.  (necause  Z  <  / 

—  k  ,  p  ~  k ,  p 


..lin  -best 

L,  +  K,  or,  Z 
k,p  k , p 


h,  <  Z  and  K  <  6 )  .  Therefore,  we 

k  ,  p  r-  k,p  k,p  = 


cannot  have  repeated  occurrence  <f  sma  1 1  changes  of  errors.  Hence,  no 


interval  of  maximal  error  can  fie  repeatedly  subdivided,  yet  its  length  not 


approach  0  . 
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Therefore,  for  all  cases,  given  arbitrarily  small  €  >  0  number  of 
iterations  p  is  finite  and  we  get  a  feasible  solution  with  criterion 
value  within  £  of  the  global  optimal  value. ] | 

4.6  Difficulties  in  Using  the  Bra  rich  and  Bound  Procedure  on  Step  Functions 
with  Decreasing  Step  Size  and  Constant  Interval  Between  Two 
Consecutive  Steps 

The  rule  of  dividing  an  interval  of  an  arc  at  an  interior  point 

(flow  level)  might,  when  applied  to  step  functions,  give  rise  to  linear 

curves  which  have  higher  values  than  the  original  cost  functions  within 

the  interval  of  interest,  as  shown  in  Figure  4.3.  So,  the  solution  using 

linear  approximation  does  not  necessari ly  give  a  lower  bound  on  the  total 

cost  and  the  procedure  is  inapplicable. 

The  above  difficulty  is  avoided  if  the  division  of  intervals  is  made 

at  the  nearest  point  of  discontinuity  (i.e.,  at  B  instead  of  D  in 

Figure  4.4).  Then  the  linear  curve  does  represent  a  lower  bound  within 

the  interval.  Each  subsequent  partition  reduces  the  error  value.  But  if 

the  restriction  (4.4c)  is  dropped,  as  it  is  for  our  algorithm,  then  for 

some  points  outside  the  interval  on  which  linear  approximation  is  based, 

the  linear  approximation  curve  lies  below  the  original  cost  function 

(shown  in  Figure  4.4).  While  working  with  interval  OB  ,  a  solution  can  be 

at  flow  level  C  and  will  have  a  positive  error  e  .  This  point  C  can 

m 

then  be  a  candidate  for  further  partition  of  interval  OB  ,  but  this  is 
impossible.  So,  for  this  method  to  be  successful,  the  restriction  of 
bounds  needs  to  be  maintained,  and  thus  one  of  the  advantages  of  our 
algorithm  is  lost.  Even  in  case  the  interval  bounds  (4.4c)  are  maintained, 
the  points  where  partition  can  be  made  are  restricted  to  points  of 


discontinuity.  So,  the  reduction  of  error  beyond  a  certain  level  will 
not  be  possible,  and  depending  on  the  structure  of  the  problem  this  can 
be  quite  high. 
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This  problem  can  be  sieved  by  using  any  arbitrary  small  €  error  by 

the  modified  method  suggested  by  J.  E.  Falk  and  R.  M.  Soland.  As  shown 

in  Figure  4.5,  where  D  is  the  point  cf  splitting  of  the  interval 

[0,u]  ,  the  cost  functions  have  open  left  and  closed  right  hand  intervals 

2 

for  each  step.  So,  it  is  a  lower  semi-continuous  function.  One  can 

c  3 

define  a  function  f  (•)  for  each  f  (•)  ,  which  is  a  convex  envelope. 

m  m 

Thus,  the  subproblems  are  convex  programs  and  of  considerable  difficulty. 

Convergence  is  guaranteed  for  any  small  €  (not  necessarily  finite 
convergence)  and  proved  in  Reference  [ F— 1 ] .  Due  to  the  unavailability  of 
a  simple  convex  programming  code,  this  procedure  is  not  tested  computa¬ 
tionally. 

4.7  Limitations  of  the  Algorithm 

The  two  major  limitations  for  this  procedure  are  as  follows: 

(a)  To  start  the  procedure,  we  need  to  specify  the  upper  bound  on 
the  flow  ii.  each  arc.  One  of  the  trivial  upper  bounds  is  the  sum  of  all 
the  flows.  This  can  be  very  large  and  if  initially  we  start  with  this 
large  upper  bound  the  error  value  is  large  and  it  will  take  a  considerable 
amount  of  partitioning  to  reduce  this  error  to  a  small  value.  Ideally, 
if  we  cor, Id  start  with  the  upper  bound  at  slightly  more  than  the  optimal 
value,  we  would  get  quick  convergence.  However,  there  is  no  theoretical 
way  to  determine  such  bounds.  At  the  same  time  if  we  assume  too  low  an 
upper  bound  the  algorithm  can  yield  a  wrong  answer,  as  will  be  evident 

^Definition:  Let  y°  e  C  .  Then  f(*)  is  lower  semi-continuous  at  y° 

if  for  every  C  >  0  there  is  a  6  >  0  such  that  if  |  jy  -  y°| |  <  6  and 

y  c  C  then  f(y)  >  f (y°)  -  «  .  f(*)  is  lower  semi-continuous  on  C  if 
it  is  lower  semi-continuous  in  each  y  c  C  . 

3 

Definition:  The  highest  convex  function  which  fits  below  f  (•)  . 

m 

A 


Cost 


)  Open  end 
1  Close  end 


Flow  <ym) 


FIGURE  4.5 
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from  an  example  discussed  below. 

A  procedure  for  generating  upper  bounds  on  arc  flows  that  wo  have 

tested  by  computation  is  as  follows:  Let  the  length  of  the  arc  m'  for 

which  the  upper  bound  is  sought  be  zero  and  let  the  length  be  f  (r  )/r 

m  s  t  st 

for  every  arc  m  /  m'  .  Here  r  is  the  minimum,  nonzero,  specified 

flow  requirement  r„  .  We  use  r^  since  if  any  flow  uses  an  arc,  the 

flow  will  be  at  least  r  and  the  smaller  the  flow  the  higher  the  cost 

st  ° 

per  unit  flow,  i.e.,  length.  Find  the  shortest  chains  between  all  pairs 

of  nodes  i  j  .  Send  flow  r..  along  the  shortest  chain  between  i  and 

ij 

j  for  all  i  j  3  r_  >  0  and  use  the  resultant  flow  on  arc  in'  as  the 
upper  bound.  While  making  arc  m*  free  and  every  other  arc  as  expensive 
as  possible  might  appear  to  give  maximal  possible  use  of  arc  m'  ,  the 
following  example  shows  that  it  may  fail  to  do  so.  It  then  might  be  hoped 
that  if  the  upper  bound  so  generat'd  fails  to  be  large  enough,  the 
resulting  solution  given  the  bound  will  be  at  the  hound  (which  can  then 
be  relaxed),  but  the  example  shows  even  this  is  false.  In  spite  of  this, 
it  seems  a  reasonable  heuristic. 


Consider  a  three  node  network  with  requirements  r^  =  3  ,  r^  =  4 
and  r23  =  ^  ar|d  cost  functions  as  shown  in  Figure  4.6. 

Upper  bounds  using  the  procedure  are  U  -2  ,  l’  -  5  ,  U  -  2  . 

J  2  1  J  c  j 

Th_  Lest  solution  respei i i ng  tiiest  hounds  is  (y  •, . v  ,  > )  ~  (1,4,1) 

with  total  cost  =  5.1.  Here  none  of  the  flow  is  at  its  upper  bounds,  so 
the  bounds  won't  lu-  relaxed.  However,  actual  minimum  cost  solution  is 
(y |  , , y |  ,.y,j)  :  ('>.<>. ’’)  with  total  «••••. t  i.h. 

Kmpiiio.il  evidence  given  l>v  :  >lvlig  the  ju  obi  inn  with  dilleieul  sets 


of  upper  bounds  shows  that  convergeiii  e  rate  for  a  giv«n  piobl 


ein  is  very 


much  dependent  on  how  small  an  upper  b  mnd  is  used  fo-  the  flows,  on  the 


Flow  (y12) - ► 


A 


FIGURE  4.6 
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arcs  of  the  network.  ft  seems  that  these  effects  are  insensitive  to 
network  structure  and  flow  pattern  (i.e.,  disimilar  networks  tend  to  show 
a  similar  ratio  of  improvement  when  bet  ter  upper  bounds  are  used) . 

(b)  The  second  limitation  of  this  method  is  the  core  size  of  the 
computer.  As  the  tree  grows,  in  a  large  problem,  we  need  to  keep  in 
store  a  large  amount  of  information  specifying  the  sets  of  intervals 

2 

corresponding  to  pendant  nodes.  Moreover,  to  take  advantage  of  only  N' 

calculations  for  the  arc  shortening,  subproblem,  we  need  the  shortest 

distances  as  well  as  information  about  Lhe  paths.  So,  for  each  node  of 

the  solut;on  tree,  three  matrices  of  size  N  x  N  (N  =  number  of  nodes 

in  the  original  problem)  need  to  be  stored.  Since  interval  bounds  arc 

integers,  we  can  pack  the  information  as  one  bit  of  information  and  use 

the  upper  half  of  a  matrix  to  store*  it,  and  use  the  lower  half  to  store 

the  distances.  So,  the  number  of  matrices  can  be  reduced  to  two  matrices 

of  size  b  *  N  .  When  k  partitions  are  made,  the  total  number  of 

matrices  is  2(k  +  1)  ,  If  nothing  has  been  blocked,  then  the  total 

? 

tiumber  of  core  locations  required  is  2N  (k  +  1  )  .  Ibis  can  be  very  large 
and  often  exceeds  the  core  size  of  a  fairly  large  machine,  even  lor  a 
network  with  only  30  to  40  nodes.  Instead  of  storing  all  the  above 
Information  in  the  core,  two  things  can  he  done: 

(i)  As  soon  as  any  set  of  matrices  is  generated,  it  is  kept  in  a 
disc  file.  In  the  testing,  an  unbuffered  binary  random  disc 
file  is  used  to  .store  tin's  informat  ion.  By  writing,  i  a<  h 
matrix  on  a  set  of  pages  and  going  directly  to  t  host  p.iges  loi 
reading  it  when  needed,  considerable  t i mo  cun  be  saved.  This 
was  attempted  with  a  36  nodi  network  using,  loin  pagt  to  writi 
each  matrix.  The  d  isadv.int  ..go  of  tin's  method  is  tii.it  o'lisiil  i- 
able  time  is  wasted  in  i  cad  tig,  and  writ  i  up.  these  matrices. 
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(ii)  Instead  of  storing  the  information  in  matrix  form  for  each 

unchecked  node  in  the  tree,  we  could  simply  keep  track  of  the 
path  and  keep  the  information  about  which  arc  is  subdivided 
at  each  node  and  at  what  point.  Then  at  any  step  all  the 
interval  information  can  be  calculated  by  tracing  the  path. 

This  requires  no  storage  on  disc  files.  The  subproblem  where 
a  particular  arc  is  shortened  cannot  use  the  modified  algorithm 
because  no  information  is  stored.  Moreover,  when  the  tree 
grows  considerably,  a  substantial  amount  of  computation  has 
to  be  done  for  generating  the  interval  information  needed.  So, 
the  saving  in  disc  file  reading  and  writing  is  offset  by 
computing  time,  and  this  method  may  not  be  very  efficient 
either. 

For  large  problems,  the  best  procedure  seems  t o  be  a  combination  of 
the  above  methods.  For  a  few  stages,  path  information  is  used  to  generate 
a  matrix  and  then  the  disc  file  is  used  to  store  the  matrices.  Thus  at 
any  point  not  too  much  calculation  is  needed  to  generate-  the  interval 
information  and  disc  reading  and  writing  is  not  too  frequent.  (The  time- 
shared  interactive  computer  program  developed  by  Stanford  Research 
Institute  named  TRKK  (Reference  [N-l])  is  good  for  this  type  of  work.) 

A  program  has  been  written  for  the  branch  and  bound  procedure  for  a 
network  of  up  to  35  nodes  and  a  number  of  problems  have  been  solved.  The 
program  and  a  brief  summary  of  data  on  solution  of  the  sample  problems 
arc  listed  in  Appendix  I!. 
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CHAPTER  5 

In  this  final  chapter,  a  post-optimization  procedure  for  problems 
with  discontinuous  cost  functions  with  decreasing  step  sizes  (Figure  1.3 
of  Chapter  1)  will  be  discussed.  Also,  a  variation  of  the  problem  in 
which  the  requirements  vary  in  different  time  intervals  will  be  formulated. 
Finally,  the  scope  of  possible,  future  research  will  be  outlined. 

3 ■ 1 _ Post-Optimization  Procedure 

Both  of  the  methods  discussed  in  Chapters  3  and  4  assume  continuity 

of  the  cost  except  at  the  left-hand  extreme  point.  However,  most  of  the 

practical  problems  encountered  do  not  have  continuous  cost  functions  but 

functions  of  the  fora  shown  in  Figure  1.3.  To  obtain  a  solution  to  this 

problem,  the  discrete  points  constituting  the  actual  cost  curve  arc* 

approximated  by  means  of  a  continuous  nondecreasing  concave  curve.  Since 

the  approximating  curve  is  noudocn-as  ing  and  concave,  it  is;  continuous; 

everywhere  except  at  the  left-hand  extreme  point  (as  shown  in  Appendix  D) . 

(Ideally,  the  continuous  nondecrcnsing  concave  cost  curve  g(y)  .should  he 

such  that,  if  f(i)  is  the  actual  cost  of  i  channels,  and  f  u  is  the 

u 

maximum  number  of  channels  possible  the  absolute  error  }  j  f  ( i )  -  g(i)| 

i  -0 

should  be  minimized  over  all  g(y)  .  However,  this  is  itself  a  difficult 
and  unsolved  problem.  So,  the  ideal  case  cannot  he  achieved,  but 
represents  an  area  of  possible  future  research.)  Then  the  method  of  either 
Chapter  3  or  Chapter  4  can  be  applied  to  obtain  the  .approximate  sulut  ion 
of  this  modified  problem  (i.e.,  with  noiuleereas  i  ng  concave  cos'.  I  fund  inn). 
Finally,  the  jol lowing  somewhat  arbitrary  local  improvement  procedure  can 
he  used  to  obtain  an  approximate  solution  of  tin  original  discrete 
variable  problem.  Some  reasons  for  choosing  this  particular  procedure  will 
he  discuss'.ed  in  the  remarks  followup,  :he  procedure. 
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Assume  there  are  n  source-sink  requirement  pairs,  and  number  these 
pairs  arbitrarily,  from  1  to  n  . 

Algorithm: 

Step  0: 

Let  i  =  1  . 

Step  1: 

Let  y  be  the  flow  in  arc  m  in  the  current  solution  to  the 
m 

problem.  Consider  the  ith  source-sink  pair.  Find  the  chain,  c;  ’ 1  it 

chain  P  ,  through  which  the  entire  flow  of  the  source-sink  pai;  i  was 

passed  in  the  solution  of  the  continuous  approximate  problem.  (Both  the 

procedures  of  Chapter  3  and  of  Chapter  4  give  rise  to  a  solutioi  with  a 

single  chain  between  each  source-sink  pair.)  Let  w  equal  the  largest 

m 

flow  value  with  cost  less  than  the  cost  of  flow  y  .  Let 

m 

t  =  f(y)-f(w).  For  each  arc  m  r  P  ,  let  7.  -  y  -  w  (see 

m  m  m  m  m  m  m  m 

Figure  5.1).  (If  the  required  flow  between  any  source-sink  pair  is 

integral  (assumed  here),  then  the  flow  level  y^  in  every  arc  is  also 

integral  (as  proved  in  Appendix  C).  Hence,  f  (y  )  is  defined.  Also, 

mm 

C  ,  defined  below  is  integral.) 
m 

Let  Z  =  Minimum  {7.  }  =  Z  ,  .  Then  reduce  the  flow  in  chain  P  lor 

,,  mm 

mcP 

source-sink  pair  i  by  amount  Z  .  This  guarantees  a  reduction  of  the 

cost  by  at  least  t  .  (possibly  more  because  there  can  be  more  than  one 

m 

arc  at  which  the  minimum  is  attained).  If  7.  is  greater  than  existing  flow 
in  the  chain  P  corresponding  to  ith  pair,  go  to  Step  3. 

Step  ? : 

For  the  resulting  flow  pattern,  determine,  for  all  m  ,  the  value, 
called  C  ,  of  the  amount  of  1  l.'W  that  can  he  sent  thiough  an  arc  m 
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without  incurring  any  further  cost  (see  Figure  5.2). 

Using  C  as  the  capacity  for  each  arc  m  ,  find  the  maximum  flow 
m 

F  between  ith  source-sink  pair  by  the  Ford  and  Fulkerson  max-flow 
labelling  algorithm  as  refined  by  Karp  and  Kdmonds  [E-l].  (The  network 
considered  here  is  undirected  and  the  labelling  algorithm  mentioned  above 
is  for  directed  networks.  So,  it  is  necessary  to  convert  the  undirected 
network  to  a  directed  one  by  replacing  each  undirected  arc  by  two  directed 
arcs  of  opposite  directions.  From  now  on,  all  chains  will  be  referred  as 
paths  because  they  are  directed.) 

If  max  flow  F  <  Z  ,  then  restore  Z  amount  of  flow  back  to  the 
path  in  Step  1. 

If  }  >  Z  ,  then  further  reduce  the  flow  in  the  path  determined  in 
Step  1  by  amount  (F  -  Z)  .  Send  amount  F  of  flow  through  the  paths 
determined  by  the  max-flow  algorithm. 

Note: 

The  C  value  for  at  least  one  of  the  arcs  in  the  patli  traced  in 
m 

Step  1  is  zero  (at  m’  where  m'  minimizes  Z  )  .  So,  the  paths  obtained 

in 

by  the  max-flow  algorithm  do  not  contain  the  previous  path. 

Step  3: 

Add  1  to  the  value  of  i  .  If  i  <  n  ,  go  to  Step  4.  If 
i  =  n  +  1  ,  let  i  =  1  and  go  to  Step  4. 

Step  4: 

If  F  <  Z  for  n  consecutive  steps,  then  stop.  If  not,  go  to 


Step  1. 
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Theorem  5.1; 

The  algorithm  described  above  terminates  in  a  finite  number  of  steps. 

Proof : 

The  finite  termination  of  the  algorithm  depends  on  the  finite  termina¬ 
tion  of  any  sequence  of  flow  reduction  steps  and  the  finite  termination  of 
the  max-flow  algorithm  in  Step  2.  Since  all  of  the  capacities  at  any 

stage  are  integral,  the  refined  procedure  of  Karp  and  Edmonds  will  yield 
a  finite  bound  of  (1  +  l°gc/c  for  each  max-flow  problem  where  c  is 

the  maximum  number  of  arcs  across  a  cut-set  and  F  is  the  maximum  flow. 

In  Step  2,  if  any  flow  change  takes  place  then  there  is  a  decrease  in 
the  total  cost  by  at  least  the  amount  of  the  smallest  discontinuity  in  any 
cost  function,  since  it  does  not  cost  anything  to  send  the  flow  through 
the  alternate  paths  found  by  the  max-flow  algorithm.  The  monotonic 
decrease  in  cost  guarantees  no  cycling.  Only  a  finite  number  of  improve¬ 
ments  is  possible  since  the  total  cost  is  bounded  below  by  zero  (i.e., 
total  cost  >  0)  .  Also,  if  no  improvement  is  obtained  n  times 
consecutively,  the  algorithm  stops.  Thus,  termination  of  the  algorithm  is 
achieved  in  a  finite  number  of  steps. ! | 

5.2  Remarks 

The  f  llowing  remarks  are  pertinent  to  the  above  procedure: 

(i)  This  is  a  procedure  to  be  used  after  a  flow  pattern  has  been 
found  by  the  method  of  either  Chapter  3  or  Chapter  4  based  on 
a  continuous  (except  at  left-hand  extreme  point)  approximation 
of  discrete  cost  data.  There  is  no  guarantee  that  we  wJ 1 1 
obtain  the  global  optimal  point  ior  the  discrete  problem.  The 
concept  of  local  optirr  il  v.  lue  is  not  applicable  here  because 
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the  variables  (flows)  are  discrete. 

(ii)  It  is  not  apparent  whether  we  could  have  obtained  more 
reduction  in  the  cost  by  considering  source-sink  pairs  in  a 
different  order. 

(iii)  For  any  set  of  source-sink  pairs,  some  flow  may  be  diverted 
from  the  initial  unique  path  to  one  or  more  paths  found  by 
solving  max-flow  problems.  Yet  Step  1  considers  only  the 
original  path  and  not  the  deviations  determined  by  previous 
max-flow  solutions.  These  other  paths  are  not  considered 
because  the  flow  in  these  paths  is  generally  small  compared  to 
the  flow  in  the  original  path.  Hence,  the  minimum  reduction 

Z  (obtained  in  Step  1)  may  be  more  than  the  flow  values  of 
these  paths.  Then,  no  cost  improvement  could  take  place. 

(iv)  In  Step  3,  if  F  <  Z  ,  then  no  improvement,  with  respect  to  a 

given  source-sink  pair,  in  the  total  cost  is  possible  using 
our  algorithm.  I t  is  possible  that  an  improvement  can  be 
made  by  using  flow  increases  greater  than  .  While  tin- 

rerouting  is  then  not  free,  it  is  possible  that  F  can  be 
made  larger  than  Z  at  a  smaller  cost  increase  than  the 
decrease  we  obtain  by  reducing  the  flow  in  I’  by  the  amount 
F  .  To  do  this,  we  have  to  consider  the  relative  sixes  of 
steps  in  the  cost  functions.  This  is  a  rather  involved 
procedure  and  was  not  attempted  here. 

(v)  The  values  of  Z  and  C  depend  on  the  width  of  the  stops 

mm 

in  the  cost  curve  (i.e.,  the  difference  between  Mow  values  at 

two  consecutive  .-.tens  in  the  cost  curve).  The  value  of  t 

1  in 

depends  on  the  height  of  the  steps  in  the  cost  function. 


There  is  not  much  improvement  possible  if  step  widths  and 
heights  are  very  small — so,  in  that  case, this  procedure  is  not 
advisable.  In  case  the  cost  functions  consist  of  a  few  big 
steps,  this  procedure  tends  to  give  reasonable  savings. 

(vi)  In  the  case  when  the  cost  function  is  similar  to  Figure  1.4, 

this  procedure  can  be  applied  with  the  redefinition  of 

and  C  as  shown  in  Figure  5.3. 
m 

5.3  Future  Research  Directions  for  This  Problem 

The  most  promising  theoretical  direction  to  pursue  in  solving  the 
general  network  synthesis  problem  with  continuous  concave  cost  functions 
seems  to  be  the  cutting  plane  method  proposed  by  Ritter.  In  this  method, 
each  local  search  procedure  should  be  a  simple  subproblem  (such  as  a 
shortest  path  problem),  and  the  generation  of  cutting  planes  to  exclude 
local  optimal  solutions  without  cutting  out  the  global  optimal  solution, 
should  be  simple  enough  to  be  handled  easily,  even  for  large  problems. 

To  make  the  procedure  of  Chapter  4  attractive,  it  is  necessary  to 
find  a  good  upper  bound  of  the  optimal  flow  in  each  arc.  Research  toward 
finding  a  reasonably  good  method  could  prove  fruitful. 

Finally,  research  seeking  better  methods  of  generating  a  good  solution 
of  the  discrete  cost  problem,  given  a  solution  of  the  continuous 
approximate  problem,may  prove  worthwhile.  But  the  problem  is  very 
difficult  and  the  procedure  will  aLmost  certainly  have  to  be  heuristLc. 

5.4  A  More  General  Research  Probl em 

A  somewhat  more  general,  and  realistic,  situation  than  discussed 
earlier  is  the  case  where  the  requirement  matrix  changes  ovi r  different 
time  periods.  The  capacity  should  be  built,  at  minimum  cost,  at  the  start 
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of  the  process  so  as  to  meet  the  requirements  for  all  succeeding  times, 
using  different  routing  during  different  time  periods.  The  programming 
formulation  is  given  as  follows: 


X^j(t)  =  A  vector  of  dimension  P  where  its  elements  x^_  (t) 

represent  the  amount  of  flow  through  the  kth  chain  between 

source-sink  pair  i  j  at  time  t  . 

A..  =  M  x  P..  path-edge  incidence  matrix,  whose  (m,k)  element  is 
ij  ij 

1  if  the  kth  chain  between  scurce-sink  pair  i  j  traverses 
edge  m  or  0  if  it  does  not. 

Y(t)  =  7  A, .X.  (t)  =  An  M-vector  where  element  y  (t)  is  the 

•  .  ij  ij  m 


total  flow  in  edge  m  at  time  t  . 

C  =  An  M-vector  where  element  is  the  total  capacity  necessary  in 

edge  m  such  that  it  meets  the  requirements  for  t  =  1,2,  ...,  T 
time  periods. 


Here  we  have  to  minimize  the  total  cost  Z  where 


M 


Z  =  l  f  (C  )  . 


m-1 


id  m 


Subject  to  C  >  y  (t)  V  t  -  1,2,  ....  T  and  for  all  arcs  m 
m  =  m 


j!  xu(t)  i  rnM  tor 

x^.(t)  >  0  for  all 
ij  • 

In  this  formulation, 


all  source-sink  i  j  and  all  t  =  1,2,  ...,  T 


t  and  l  j  . 


the  number  of  variables  has  increased  enormously 


An  efficient  solution  procedure  is  not  in  sight  and  can  only  be  found  if 
the  problem  of  this  thesis  is  solved  first. 
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However,  an  even  more  complex  problem  faces  the  telephone  company. 
The  requirement  matrix  changes  over  succeeding  years  and  can  be  forecast 


at  the  beginning  of  the  planning  period.  The  cost  of  rearrangements 
(i.e.,  dismantling  the  facilities  between  certain  nodes  and  putting  them 
between  certain  other  nodes  in  different  periods)  of  facilities  is  also 
given.  An  initial  network  has  to  be  synthesized  which  will  meet  the  initial 
capacity  requirements.  Additions  and  rearrangements  in  succeeding  years 
must  be  planned.  There  may  also  be  yearly  budgetary  constraints.  All  of 
these  must  be  solved  in  such  a  way  that  discounted  total  cost  over  several 
years  is  minimized.  This  problem  is  much  more  complex  than  the  one 
attempted  in  this  thesis  and  a  heuristic  solution  procedure  can  be  a 
challenging  extension  of  this  work. 
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APPENDIX  A 


AAS128  22/27/72 


5C  LANGUAGE:  FORTRAN  FOUR,  GE  635  (MARK  TWO)  COMPUTER 
I2C  132  K  MACHINE,  1  MICRO-SECOND  MEMORY  CYCLE 
15C  SIMILAR  TO  13M  36.0-65  OR  CDC  6422. 

22C  PROGRAM  FOR  A  LOCAL  SEARCH  ALGORITHM  OF  SOLVING  A 
25 C  MULTI-COMMODITY  FLOW  PROBLEM  UNDER  CONCAVE  COST. 

32  C  GIVEN  FLOW  REQUIREMENT  MATRIX  R,  COST  MATRIX  K 
35C  COST  OF  F ( I  , J )  FLOW  ON  ARC  ( 1 , J ) =K C I  , J  ) *F ( I , J ) **? . 
43  C  N- NUMBER  OF  NODES  IN  THE  NETWORK. 

45C  INPUT  DATA  FILES  SHOULD  SUPPLY  VALUES  OF  R  *  K. 

50  FILENAME  2DATA 

55  INTEGER  R ( 34 ,3 4 )  ,CM (34  ,3 4 )  ,C (3 4 , 34 ) , F (35 , 3 5 )  ,S 
60  INTEGER  Y 

65  REAL  K (34 ,34)  ,0(34,34)  ,C0(34, 34  ) 

7 2  10  FORMAT  (  5X  ,  F 0 , 2 , 1  5  ) 

75  20  FORMAT  (5X.F8.3) 

SO  32  FORMAT  (13  ,4X,  13 ,14) 

85  4  0  FORMAT  ( 3X  ,13) 

9(1  52  FORMAT  (213,14) 

95C  SUPPLY  THE  INFORMATION  ABOUT  N,?  ,  INPUT  FILE, 

1 2 PC  VALUE  OF  COEFFICIENT  OF  CONVEX  COMBINATION  ’3 
125C  AND  INITIAL  NUMBER  OF  SAMPLE  NO. 

1  ID  PHI  NT  ,"  N ,  Z  ,  Z  D  AT  A  ,3  ,  NO" 

115  I  NPUT  ,  N  ,Z  ,Z  D  AT  A  ,B  ,  NO 
122  NUN-1 
125  DO  1  OK  1  =  1,  N1 
132  I  1=1  +  1 
135  DO  1  »?l  J  =  I  1  ,N 

143  122  READ  (ZDATA,13)  K  ( I  ,J  )  ,R  ( I  ,  J  ) 

145  PRINT,  I ,  J  ,  K  ( I  ,  J  )  ,R  ( I  ,J  ) 

15fC  STARTING  FROM  A  RANDOM  POINT 
155  S - 2  ;  I T R A C  =  0 
162  152  S=S+1 

165  DO  1  7,1  1  =  1  ,N1  ;  II  =1  +  1 
172  DO  1  72  J=I  1  ,  N 
175  IF(K(I,J).LT.  12  22)  GO  TO  167 
182  D ( I , J ) = 1232  ;  GO  TO  17/ 

185  167  Y=RRANO(  I  .3)  ;X=RND(1.2) 

U3  F  ( I ,  J  )  =  (  122  *v+  1  2*Y  )  /P 
195  PRINT  5f,I,J,F(I  ,J) 

222  D(I,J)  =  (7.*XU,J))/CF(I,J)**(  l-Z)) 

225  1 72  D(J , I )  =  D( I  ,J) 

212C  CREATION  OF  INITIAL  CH(I,J)  MATRIX 
215  262  DO  2  73  I  =  1  ,  N 

223  DO  2  73  J  =  1  ,  N 
225  272  CH  ( I  ,J)=J 
232  DO  32  3  U  1  ,  N 
235  332  0(1,1)  =  2 

2420  STARTING  OF  ITERATIVE  PROCESS. 

245  I T  R  A  =  3 ; N  C  3  V  =  8 ; NC  0  =  2 
25i*C  CREATION  OF  C(I,J)  MATRIX 
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29?  IF (0(1  ,M)  .GE.  1000)  GO  TO  500 
295  DO  500  J  =  I  1  ,  N 
300  A  =  DCI  ,  WD+DOI.J ) 

305  IF  C!)(I,J).LE.A)  30  TO  500 
31  0  Dd,J)=D(I  ,X)+DC%J) 

315  D(J,I)=D(I,J) 

320  C  (I  ,J)  =  C(I  ,w|) 

325  C(J,  I  )=C(J  fM) 

330  5 flfl  COMTINUE 

33 5 C  CHECKING  FOR  CONVERGENCE 

340  ICKE=0 

34  5  DO  520  I  :  1  ,  N 

350  00  520  J  =  1  ,  N 

355  IF  (CM  ( I ,  J  )  .E 0.  C(I,J)>  GO  TO  510 
350  CH(I  ,J)=C(1  ,J> 

365  GO  TO  520 
370  510  ICHE  =  ICHE+1 
375  520  CONTINUE 
380  NCH  E  = N*  N 

385  IFdCHE.EQ.  NCHE)  NC0  =  2 

390  I F ( I  CHE  .£Q . NCHE  .AND.  S.E3.1)  NC0V=2 

395  IFdCHE.EQ.  NCHE  .AND.  C  OG  TT .  GT.  C  OCT )  \'C0V  =  2 

40 (PC  UPDATING  THE  FLOW  MATRIX 

405  DO  550  Is  1,  VI 

410  11=1+1 

415  DO  550  J  =  I  1  ,  N 

420  550  F(  I  ,J)  =  0. 

425  IFCNCOV  .  EQ.  O)  GO  TO  555 
430  PR  I  NT,"F  R01.BT  ,  TO,  FI  OW" 

435  555  00  6P?  I=1,NI 

440  11=1+1 

44  5  DO  500  J  =  I  1  ,  N 

450  IF(R(I,J).E3.0)  GO  TO  60? 

455  IF  CNCOV.  NE.  P)PRINT  30 , 1  ,J  ,R  ( I  ,J  ) 

450  X  =R ( I  ,  J  ) 

455  Y  =  I 

47 a  550  L  =  C( I , J ) 

475  IF  (I . GT.L )  GO  T  1  570 
480  F(I ,L)= (I fL)+X 
485  GO  TO  58/ 

490  570  F(LtI)=F(L,I )+X 

495  580  IF  (L.EQ.J)  GO  TO  590 

500  IFCNCOV.  NE.fl)  PR  I  N  T  40  ,L 
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I=L  //‘W? 

GO  TO  5 SB 
SO"  I  =  Y 
6/B  CONTINUE 

IF  (NCO  .FQ.  2)  GO  TO  75/ 

IF  (NCOV  .  EQ.  i)  ■JCOV  =  H 
CALCULATION'  OF  NEW  DISTANCE  MATRIX 


5/5 
01/ 

51  5 

52  Z 
525 
530 
535C 
54/  COST:/ 


tt 


545  DO  783  1  =  1,  01 

55/  11=1+1 

555  DO  7/3  J  :  I  1  ,  N 

56/  CO  (I  ,  J  )=  (K  (I  ,J  )  )*(  C  F  ( I  ,J)  >**Z> 

56  5  COST  =C  OS  r+CO<I  ,J  ) 

57/  IF  (F ( I ,J ) . EO.B)  F( I ,J )  =  1 

575  D(  I ,J )  =  C (Z*K ( I ,J ) ) /(F( I ,J) +  * (1 -Z ) ) )*9+( 1 -3 ) *  CD C I  , J  )  ) 
58B  D ( J  ,  I  )  -OCX  ,J) 

5P5  7BB  CONTINUE 

59B  PRINT, "TOTAL  COST:", COST 

505  I  TR  A  :  I  TR  A+  1 

6/3  GO  TO  4 rr 

6/5  75/  PRI NT  ."CONVERGENCE  AT  ITERATION  NO  =",ITRA 
61 B  I  TRAC : I TR  AC+ 1  TR  A 
615  IF  (S.GT.  I  )  GO  TO  78/ 

62/  C  OS  T  T  : C  OG  T  :  COSTU:COST 
62  5  GO  TO  1  r 

63 B  78/  IF  (C Oc TT  .GE.  COST)  COSTT:COST 
635  IF  (  COST'J  .LT.  COST)  ''.OS TO: COST 
640  IF  (S.LT.MD)  GO  TO  15/ 

645  ERR:  (C OSTU -COSTT  )  /  CG  -  1  ) 

653  IF  CS.NE.NO)  GO  TO 


655  CM: I  TRAC/ S 

66  .*  r  t ,  I  N  7  ,  "  ER  R :  "  ,  f.R  R  ,  "C  M  :  "  ,  C  '1 

66 5C  SUPPLY  TWF.  COEFFICIENT  OF  THE  LOS''  FUNCTION  AT) 
67/C  C  .IMPUTATION  COST  FUNCTION. 

67t.  PRINT,  "  r. ,  C  N  " 

60/  J  npu:  ,  V. ,  0  N 

60  5  0  5/  ;  F  ( (  •>  E:<  R  )  .  L  T .  (C  N*  CM  ) )  So  TO  I  r"? 

60/  NM:  NO-*  ^ 

600  IF  (  '■  .ST.  NM)  SO  T  o  1/0/ 


7/ 1*  S  )  TO  153 


fr  5  1  3  0  ?  PR  IN  r  ,  C  :)S  T  r  ,  s  ,  F.KR 
71/  '.T  OP  ;  Z N : > 
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5C  LANGUAGE:  FORTRAN  FOUR,  GE  635  (MARK  TWO)  COMPUTER 
1BC  132  K  MACHINE,  1  MICRO-SECOND  MEMORY  CYCLE 
15C  SIMILAR  TO  ISM  368-65  OR  CDC  6488. 

2CC  PROGRAM  FOR  A  GL03AL  SEARCH  (BRANCH-BOUND)  ALGORITHM 
25C  OF  SOLVING  A  MULT  I - COMMOD I T Y  FLOW  PROBLEM  UNDER 
3fC  CONCAVE  COST.  GIVEN  FLOW  REQUIREMENT  MATRIX  R,COST 
35C  MATRIX  KS,  N=NUM3ER  OF  NODES  IN  THE  NETWORK,  COST  OF 
ABC  FY(I,J)  FLOW  IN  ARC  ( I ,J) =KS( J, I)+KS( I , J) *FY( I , J) *** 
45C  INPUT  DATA  FILE  SHOULD  SUPPLY  VALUES  OF  R  A  KS 
COMMON  DY,DX,FY,R,CY,CX,N,NE,K,L 
50  REAL  DY(35,35),DK(35,35),ZA( 188), ZL( 188 ),KS(35,35) 

68  REAL  D ( 1268) 

65  INTEGER  S  (  1  88  )  ,  KE  (  1  88  )  ,LE  (  1  88  )  ,CY  ( 3  5 , 35  )  ,CX  ( 3  5 , 35  ) ,  RM 
7 8  INTEGER  CU268) 

75  INTEGER  R  ( 3  5 , 3  5  )  ,  FY  (  3  5 , 35  )  ,L3  ( 3  5 , 35  )  ,F  (  !  88 ) 

88  FILENAME  ZD  AT  A 
85  If  FORMAT  ( 5X , F9 .2 , 1  5  ,  F8  .2 ) 

98  28  FORMAT  (13,  AX, 13, 14) 

95  3 8  FORMAT  (3X.I3) 

108  32  FORMAT  (I  7,  FI  8. 2) 

185  3A  FORMAT  ( I  7 ,2F l « .2) 
llfC  SUPPLY  THE  VALUES  OF  N  ,7. ,  I NPUTF1LE  ,  NUMBER  OF  PAGES 
115C  FOR  EACH  DATA  BL  OCK ( = NUMBER  OF  V AR I ASLES /3 1 5  ) 

128  PR T  NT,  "N ,Z , Z  DAT  A , NP" 

125  INPUT,  N  ,Z  ,Z  DAT  A  ,  NP 
138  M  1  =  M-  I  •,  RM:  18*88 
135  DO  18  8  I  =  1 , N 1 ;  I  I  =  1+  1 
148  DO  188  J  =  I 1  ,N 

145  READ  (ZD  AT  A  ,  18)  KS  ( I  ,J  )  ,R  (  I  ,  J  )  ,K  S  ( J  ,  I  ) 

158  Iff  IF  <  (R  (  1  v.)  >  .ST.  8)  .  AND.  (R  ( 1  ,J)  .LE.RM)  )RM=R(  I  ,J  ) 

155  DO  158  I  :  1  ,  N1  ;  11=1  +  1 
168  DO  1  58  J  =  I  1  ,  N 

165  IF  (KS(I,J).LT.  1  8  8  8)  GO  TO  128 
178  DY( I ,J>= 1888 
175  GO  TO  158 

1  SB  128  OY(I,J)  =  ((KS(J,I)+KS(I,J)*(RM**Z))/RM) 

185  158  D Y ( J , I ) =  DY ( I  ,  J ) 

198  DO  168  1  =  1,  N 
195  DYd  ,1  )  =  K 
288  168  CX(I,  I) =1 
28  5  ME  zf 

218  CALL  SHORTEST 

21 5C  OBTAINING  BOUNDS  ON  FLOW  ON  EACH  ARC. 

228  DO  388  K  =  I  ,  Ml  ;  K  1  =K+  1 
22  5  DO  388  L=K  I  ,  N 

23?  IF  (KS(K,L).GE.  1888)  GO  TO  298 
235  DX  (K  ,L  )  =  8 
248  CALL  ARCHANGF. 

245  R(L  ,K  )  =Y(K  ,L  ) 

258  DYCK  ,L)  dKS(L  ,K  )  +  (KS(KfL)  )*((«KL,K)  )**Z ) ) /R (L  ,X  ) 


112 


A  A  S3??  1  r,2/r.?/ 72 


255  G'J  TO 

20/  ?.3?  0  Y  ( K  ,  L )  - 1  “ » 

255  3/7  0  Y ( L , K) r  0  Y ( K , L ) 

27rC  START  I  Vo  Of  3R4VCH  -t  tJOMVD  PR 0CE9URE. 

275C  CREATIVfi  A  R  A  VO  01  31  VARY  FILE  TO  STORE  THE  OAT  AS . 

2 P.r  OPEVFILE  SOATA?**,":?",  1  t252«.V> 

285  CALL  IF  3REAK ( >753) 

293  VE=l  ;  -IF =  I 
295  V'JIrl 
31’ 7,  CALL  SHORTEST 
3?5  00  33-*  I  r  1  , 1  rn 
31 r  3  3/  S<I):? 

3 !  5  S  (  I  )  r  1 ;  1 1  =  I 
32 /  90  3  5/  I -I  ,  VI:  11  =  1+1 
325  00  35/  Jr  1 1  , V 
332  L2  ( I  ,J)=r?  (.1,1) 

33  .5  35°  L  3  ( J  ,  I  )  =  7 
3 A ?.  SO  TO  5;’ 3 

34  5  4/1  X  =  K  F  (  V  F  )  •  L  =LF  C  VF  ) 

352  VI)  Or  VJ  3+1 

355  IF  (UF.EQ,Y£>  30  TO  42-' 

3C  7  1  P  :?.+  (  Yr  -  l  )  *  \IP+  1 
3Go  CALL  UREA)  C  S  0  A  T  A  ,  ) ,  I  P ,  JP  ) 

3  7/  I  2  r I P  +  V  p 

37  0  CALL  2EEA9  (  S  )  AT  A  , C  ,  I  p  ,  V'1 5 

383  VO  =  V 

325  0  0  -115  V  1  ,  '.'U 

35  7  J  r  I  V  T  (  (1-1  )/V)  +  l 
395  Irl-  V*  ( J  -  |  ) 

4?r  CYU  ,J)  =  CM) 

A!'  5  IF  (I.LE.J)  SO  TO  A  1  J* 

Air  LKJ,I)rIvrC)(1)) 

415  L-3U  ,  J)  =  (  )(  1)  -L  3  ( J  ,  I)  )*!'•*/ 

42/  SO  TO  4  1  3 

425  112  9  v( i , J ) r  )(1) 

43/  OY(J,  I  )  r')V(  ;  ,J) 

43  5  415  C'miVUi 

14/  42/  )/(/,!.)  :KS(K  ,L  )♦  CL  .UY  ,L.)  >*Z-F(  VF)**7  )/  (L  UX  ,L)-F  ("D  ) 
445  S  ALL  APvS’AUOf 
45/  L  > L K  r L 3 ( L  ,  K  ) 

455  L  U! ,K  ):F('.V1 
45/  00  43?  I r ! , n 
45  5  4  3/  IF  (5(  15  SO  TO  14? 

4  7/  4  1/  V E r  J  ;  SC.'Drl 
475  JO  TO  55/ 

48 r  5 7  ?  Or  r  OF 
4S5  VUI r  \"JVH  1 

40,’  5?5  0  0  5  4/  1  :  1  ,M  ;  I  I  r]  +  1 
405  0  0  54  7  Jr  II  ,  V 

5/1*  IF  (KS(I  ,J)  .0:1.  ir.f-3)  5.0  TO  53/ 
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5?  5  IF(LJ(J,  I)  .'ZO.OGD  TO  5  I  / 

51  r  OY(  I  ,J)  =KG(  1  , J  )  «(L0(  1  ,2  )  +  *3-L3(  J  ,  1)  +*?)  /  (L  HI  ,  J)-L°.  ( 
515  30  TO  54>". 

52!’  5  1?  0Y(  I ,J5  -  (KO(J  ,  I)  4K?(  1  ,.J  >*<L3U  ,J  >**!)  )  /ri(l  ,.l  ) 
525  GO  TO  5 A? 

53?  5  3?  OYU  ,J)-  ! 

53  5  5-1?  OYdJ ,  I )  -  OY  ( I ,  J  ) 

51?  IF(l,-i£.  EO.  I  )  3  )  TO  510 
515  IF  (L3LX .70.2)  GO  TO  5  13 

550  OY(K  ,L  ,L)  H  (F(  'M)**Z  )-  CL3LK  +  *7  )  )  /  (FC;-)-’  3L:/  ) 

555  GO  TO  5  16 

56 f  513  OY(K  ,L)  -  (KG(L  ,K)+KG(K  ,!_)*(r  (  ) )  /F  (•;•-.) 

56  5  515  JYO.  ,X  >-OY(K  ,L) 

5  70  LJlCK  ,L>  =  <  VF) 

5  75  L.i(L  ,K  )  :L5LK 
50/  5  10  CALL  03  0 
505  i  h  r  r 
50/  3  5/  0  0  55  3  I 
50  5  0  0  5  35  J  -  I  1  , 

50?  1-(J-  1  )•*';+ 1 

60  5  0(0)  ;  >Y(  I  ,J) 

61  r  if  (Fyc  ,J) .LF.r;  (j  ,1 ) )  00  to  555 

515  IF  (FYU  ,J).L7.L  HI  ,  J  )  )  0)  y.j  353 
52'-‘  L  i(  I  ,J  ):7Av:-(r'Y(  I  ,.J  )  ,2*  ,1  ,  1  )  1 

S3 5  11  :  (I-  1  ><  *,:+  J 

53?  0(11  )  rL .1(1  ,J>  +  CL.K.J  ,  I  )  /  l  '?•’  ) 

635  1  :  1 

51 0  355  CO  or; 

sip  1F(1.{7. 70.  1  )  GO  f  0  5  ■'  3 

650  00  550  1:1,1 

65  5  J  O  55 ?  J  -  1  ,  5 

66?  (J- 1  )*■'+ 1 

665  5G?  C  (  DiCYU  ,J  ) 

5  7?  I  P;2*  (  OF-  1  )  *‘ip+  I 

675  CALL  0001  fF  (  1  1  A  I  A , 0  ,  I J , OF ) 

SO?  I  P  ; I  P  +  ‘J p 

603  CALL  GV.arF  (CO  A  F  A  ,r  ,  1  P  ,  \'P  ) 

GO/C  COOT  CALdLATIOVO 
60  5  6/--  ?  L  (  V  ‘  l  7  A( 

7  /  <  0  0  6  5<*  I  -  1  ,  V 1  :  I  I  z\  +  1 
7 f  5  OO  55?  J :  I  1,0 

71/  1FU.H1.J)  .OF.  1  /?  r )  00  TO  65-' 

71  5  7  A  AH  HJ  ,  I  )  1  K~(  I  ,J  )  ‘  (FV(  I  ,  I  )  »  « "  ) 

72/  I  r  (L‘i(  J  » ! )  .  70.  r )  GO  T  1  6  'r 

72  5  ZLL-  *  Y  (  ',..')*(  FY(I  ,2  (J  ,  1  )  )  '  7  0  f,i ,  n  -tX  5  (  J  ,.l)  H  ’  H.l 

73/  GO  TO  61? 

733  03?  7LL:  )Y(  I  , J  )  *FY(  I  ) 

71?  61/  1  A(  O?;)  z~  A(  G  7)  +-  AA 
713  ?L  (  V£)  : 7  L  (  +ZL.. 

75/  0X(J,I):7AA--LL 


J  ,  1)  ) 


,  n  *  *- ) 
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755  IF  <ERR.GE.DX(J,I>>  GO  TO  S5f 
760  ERR: DX  (J  ,  I ) :  K E(  NE)  =  I:  LE (  NE)  =  J 
765  FCM£)  =  FYU,J> 

77f  550  CONTINUE 
775  IFCII . N E. 1)  GO  TO  670 
780  PRINT, ERR.ZAC" E) ,ZL(NE> 

785  PRI  NT  ."VALUE  OF  ERROR  PERMITTED  EP” 

790  INPUT, EP 
795  Z9=Z  A (  1 ) 

800  S7e  IF  (Z  A(  NE) . GT  ,Z3)  GO  TO  630 
805  Z8:ZA(NE);N3: ME  :ZC=Z3-EP 
81 0  PRINT  32 , N3,Z3 
31 5  680  IF(NE.EQ.Il)  11:11+1 
820  PRINT  34 , NE,ZA(N£) ,ZL(NE) 

82  5  IF(NE.NE.NF)  GO  TO  50.0 
830  DO  70.1  I:  1  ,  II 
835  IF  (5(1)  .£0.0)  30  TO  74/5 
840  IF  (ZL(I).GE.ZC)  S(I)=0 
345  70/  CONTINUE 
850C  CHECKING  FOR  TERIl NATION 
355  CHEC:f 
860  DO  710  1=1,11 

865  710  I F ( 3 ( I ) .EG. 7)  CHEC:CHSC+1 
370  IF  (CHEC. EQ. II )  GO  T 0  75 0 
875  Z LI  =ZL(  I  ) 

880  00  720  1:1,11 
885  I F (5 ( I ) . EG. 0 )  GO  TO  720 
890  I F (ZL  C I ) .GT.7LI )  GO  TO  720 
895  Z  L I  =Z  L  ( I ) ;  NF  :  I 
900  72 0  CONTINUE 

9U  750  PRINT, "LEAST  COST  5  OL'JT  I  ON:  "  ,ZB  ,  "ERROR  PERCHED:"  ,  EP 
915  PRINT, "NUMBER  OF  EVALUATION  :  ",  NU>) 

92  0  IP :2* ( N3-  1  ) *  NP+  NP+  1 

925  CALL  URF.AD  (S  0  A  T  A  ,  9 , 1  P  ,  NP  ) 

930  PRINT, "FROI.^T  , TO, FLOW" 

935  DO  755  >18  1  ,  N'J 
94f  J:INT(  (>1-1  >/N)  +  l 
945  I:>1,  N*  (J-  1  ) 

950  755  CY(  I  ,J)  :C(>1) 

955  DO  7 80  I:  I  ,  Mi  ;  I  1  :I+  1 
960  DO  780  J:I  I  ,  N 
965  IF (R(l , J 0 . EG . 0 )  GO  TO  784 
970  PRINT  20,  I  , J  ,R  ( I  ,J  ) 

975  NYd 

98f  760  L  0:  CY  (  N Y ,  J  ) 

985  IF  (LO, EQ. J)  GO  TO  7g0 
950  PRINT  30,  LO  ;  NY=LO 
995  GO  fO  760 
1000  780  CONTINUE 
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1005  STOP  ; END 
1010  SUBROUTINE  SHORTEST 
1015  COWON  dy,dx,fy,r,cy,cx,n,ne,k,l 
1020  REAL  DY(35,35),DX<35,35> 

1025  INTEGER  M ,  NE  ,CY  ( 3  5  ,3  5  )  tCX  C  35  ,3  5 )  ,FY  <  35 ,3  5  )  ,R  ( 3  5 , 35  ) 
1030  Ml  =N-1 
1035  DO  850  1  =  l  ,  N 
1040  DO  850  J  =  1  ,  N 
1045  850  CY<I,J)=J 
1050  DO  862  y\z  1  ,  M 
1055  DO  862  1  =  1  ,  N1 ;  II  =1  +  1 
1060  I F ( DY  C I , X) . EQ. 1 00  0)  GO  TO  862 
1065  DO  860  J=I 1 ,N 

1070  IF(DY(I,J).LE.COYU,M)+DY(M,J)))  GC  TO  850 
1075  DYCI,J)=DYCI,M)  +  DY<!*1»J) 

1080  DY(J,I)=DY(I,J> 

1085  CYCI,J)=CYCI,M);  CY( J , I > =CYC J , 1) 

1090  860  CONTINUE 

1095  862  CONTINUE 

lire  IF(NE.EQ.O)  RETURN 

1105  DO  870  1  =  1 t  N 1 ; J 1 = I + 1 

1110  DO  870  J  =  I 1  ,  N 

1115  870  FYCI,J)=0 

1120  DO  900  1  =  1  ,  N1  ;  I  1  =  1+  l 

1  125  DO  900  J=I  1  ,N 

1130  1F(R  ( I  ,J  )  ,  EQ.O)  GO  TO  900 

1135  NX=RU,J>  ;  N  Y=  I 

1140  880  LO=CYCNY,J) 

1145  IF(NY.GT.LO)  GO  TO  885 
1150  FY(NY,LO)=Y(NY,LO)  +  NX 
1155  GO  TO  890 

1160  885  FYCLO,  NY)  =  FY(LO,  NY)  +  NX 

1165  890  IF(LO.EQ.J)  GO  TO  900 

1170  NY=LO 

1175  GO  TO  880 

1180  900  CONTINUE 

1185  RDTURN 

1190  END 

1195  SUBROUTINE  AR CHANGE 

1200  COWON  DY,DX,FY,R  ,CY,CX,N,NE(K,L 

1205  REAL  DYC35.35) ,DX(35,35) 

1210  INTEGER  N  ,K  ,  L  ,C  Y  (  3  5  ,3  5 )  ,CX  C  35 ,3  5 )  ,FY(  3  5  ,3  5  )  ,R  C  3  5 , 35  ) 
1215  NUN- 1 

1220  DO  920  1  =  1 ,  N1  ;I  I  =1+  I 

122  5  DO  920  J  =  1  1  ,  N 

1230  AD=DY(I,X)+DX(K,L)+DY(L,J> 

1235  IF  <DY(I ,  J)  .LE.  AD)  GO  TO  910 

1240  DXCIfJ)=AD;CXCI,J)=CY<I ,K);CX(J,I)=CY(J,L) 

1245  IFCI.DQ.K)  CX(I,J)=L 
1250  IF(J.EQ.L)  CXCJ,I)=K 
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1255 
126(5 
1265 
127  0 
1275 
128E 
1285 
1290 
1295 
138  8 

130  5 

131  8 
1315 
1328 
1325 
1338 
133  5 
134(5 
1345 
135  0 
1355 
1368 
1365 


GO  TO  92(5 

918  AD  =  Dy(I  ,L  )+DX  (K,L)-i-DYCKfJ) 

IF(DY( I ,J) ,LE. AD)  80  To  9J5 

DX Cl ,J ) :AD; CX(I  ,J  )  =CYC 1  ,L) ;CX(J,  I  ) rCY( J  ,K ) 

IF(I.EQ.L)  CXCI,J)=K 

IF(J.EQ.K)  CX(J,I)=L 

GO  TO  928 


1 1  =  H  1 


;  1 1 =1+  1 


928  continue 

00  938  Is  I  #R|1 

DO  938  Jsl  1  ,  N 
938  FYCI,J)  =  8 
DO  988  Is  1 
DO  988  J=I  1  ,  N 
IF  (R  ( 1  ,J  )  .  EQ.8)  GO  TO  98 
NX=R(lfJ)j  MY=I 
948  LO-CX(NY,J) 
IFCNY.GT.LO)  GO  TO  958 
FYC  MY,L  0)  =  FYC  NY,LO)  +  NX 

GO  TO  968 

958  FYCLO,  NY)=Y(L0,,NY)  <  MX 
960  IF  CLO.EQ.J  )  GO  TO 
NY  =  LO 


1370  GO  TO  94  4 

1375  980  DYCI,J):0XCI  ,J  ) 

1380  DO  925  1  :  1  ,  N 
1385  DO  985  Js  1  ,  N 
1390  CYCI,J)  =  CX(I,J) 

1395  9«5  IF  CI.LI.J)  DY C J , I) = DY C I fJ ) 
1408  RETURN 
1405  END 


THE  COMPUTATIONAL  RESULT  OF  A  NUMBER  OF  PROBLEMS 
IS  LISTED  INTI  IE  TABLE  OF  NEXT  PAGE. 


TABLE  OF  COMPUTATIONAL  RESULTS  FOR  THE  GLOBAL  SEARCH  PROCEDURE 


CRU  is  approximately  1  second  of  computing  tine. 
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APPENDIX  C 

The  solution  procedures  considered  for  this  problem  always  use 
continuous  functions  and  continuous  variables.  However,  except  for  the 
solution  method,  this  assumption  of  continuity  of  the  cost  function  is 
not  needed.  The  values  of  the  cost  function  at  integral  points  are 
enough  data  to  get  the  total  cost.  If  the  requirements  (r  )  are 
integers,  then  the  following  observation  shows  that  all  the  extreme  points 
(one  of  which  is  found  at  each  iteration  of  our  procedures)  of  the  convex 
polyhedron  are  a1cj  integral. 

The  problem  is: 


M 

Minimize  Z  =  7  f  (y  )  , 

,  m  m 
m=l 


k  k 

subject  to  x. .  >  r.,  V  (i  j  nairs)  and  x..  >  0  ,  where 
J  ij  =  lj  '  ij  = 

l  A. .X.  .  =  Y  . 
all  i,j3i>j  1J  1J 


This  is  equivalent  to: 


Minimize  Z 


(X)  =  l  f  (  l  am.X..\ 
,  ml  .  . L .  .  i  j 
m=l  J  / 


subiect  to  DX  >  P.  ,  x.  .  >  0  ,  where  a..  is  the  mth  column  of  the 
J  =  ij  = 


matrix  A, .  and 
ij 


D  is  a  f  \  P.  .1  *  n  matrix  (n  =  the  number  of 

V1 » J  1J/ 


source-sink  pairs)  with  P  l's  in  each  row  and  R  =  {r„}  »  an 
n-vector.  Now, 
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111 


11111 


is  a  totally  unimodular  matrix. 

The  theorem  of  Hoffman  and  Kruskal  [Page  125,  Ref.  H— 3 ]  proves  that 
if  the  vector  R  is  integral  then  all  the  extreme  points  of  the 
polyhedron  defined  by  the  constraints  are  integral. 


^Definition:  A  matrix  is  said  to  be  totally  unimodular  if  and  only  if 
every  subdeterminant  of  D  Is  e<iunl  to  +1  ,  -1  ,  or  0  . 
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APPENDIX  D 

The  property  of  concave  nondecreasing  functions  established  below  is 
used  in  Chapter  3. 

Theorem: 

Let  f(y)  be:  (i)  a  concave  function  in  y  e  [0,a]  ,  (ii)  f (y) 
(i.e.,  it  is  bounded  below)  and  (iii)  nondecreasing  in  y  ,  then 

(a)  f(y)  is  continuous  (except  perhaps  at  0)  .  If  it  is  not 
nondecreasing,  then  it  is  not  necessarily  continuous  at  point 
a  . 

(b)  Both  the  left  hand  and  the  right  hand  derivatives  D  f  and 

■f-  — 

D  f  exist  at  all  points  (except  at  0  where  D  f  is  not 

defined  and  D+f  may  not  be  finite;  and  at  a  where  D+f  is 

not  defined). 

(c)  The  following  inequalities  are  true: 

(i)  D~f(y1)  >  D  f (y 2 )  v  y2  >  yl 

(ii)  D+f(y1)  >  D+f(y2)  V  y2  >  y1 

(iii)  D~ f (y )  >  D+f(y)  V  y  E  (0,a)  . 

Proof : 

This  proof  has  been  done  by  Hardy,  Polya  and  Littlewood  [Page  91, 
Ref.  H—  2 ]  for  any  convex  function  defined  in  the  open  interval  (H,K) 

t 

which  is  bounded  above  in  some  interval  i  interior  to  (H,K)  If 


The  inequalities  (i),  (ii)  and  (iii;  are  the  reverse  of  the  inequalities 
in  the  theorem  of  Hardy,  Polya  a>’d  Littlewood  because  they  are  for  convex 
functions. 


AH 
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f(y)  is  a  concave  function,  -f(y)  is  convex.  Condition  (ii)  guarantees 
a  lower  bound  for  the  concave  function  which  is  equivalent  to  the  upper 
bound  of  the  convex  function  -f(y)  .  Here  (H,K)  =  (0,a)  .  So,  the 
proof  is  exactly  similar.  Condition  (iii)  guarantees  that  at  a  there 
cannot  be  a  discontinuity  because  then  the  value  of  f  would  have  to 
decrease. 

Thus,  this  theorem  is  valid  for  a  piecewise  linear  concave  function 
(a  continuous  version  of  the  function  as  shown  in  Figure  1.5). 
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